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S UMMARY 

This document provides logical flow and guidelines for the construction 
of a low thrust orbit determination computer program. The program, 
tentatively called FRACAS (Filter Response Analysis for Continuously 
Accelerating Spacecraft) , is capable of generating a reference low thrust 
trajectory; performing a linear covariance analysis of guidance and 
navigation processes, and analyzing trajectory non-linearities in Monte 
Carlo fashion. The choice of trajectory, guidance and navigation 
models has been made after extensive literature surveys and investigation 
of previous software. A key part of program design relied upon experience 
gained in developing and using Martin Marietta Aerospace programs : 

TOPSEP (Targeting/Optimization for Solar Electric Propulsion) , GODSEP 
(Guidance and Orbit Determination for SEP) and SIMSEP (Simulation of SEP). 




1. INTRODUCTION 


A major requirement for spacecraft systems design is an effective 
analysis of performance errors and their impact on mission success. 

This requirement is especially necessary for low thrust missions where 
thrust errors dominate all other error sources. Fast, accurate parametric 
error analyses can only be performed by a computer program which is efficiently 
constructed, easy to use, flexible, and contains modeling of all pertinent 
spacecraft and environmental processes. The FRACAS (Filter Response 
Analysis for Continuously Accelerating Spacecraft) program is designed 
to meet these characteristics. It is intended to provide rapid evaluation 
of guidance, navigation and performance requirements to the degree necessary 
for spacecraft and mission design. 

This document describes the structure of FRACAS. The three basic 
program modes (trajectory generation, error analysis, simulation) are 
integrated in a master program which selects appropriate routines and 
performs the necessary executive control. The total primary and secondary 
overlay structure will require less than 70,000 octal words of a CDC 
6000 series computer. Descriptions of the overall logic (macrologic) and 
of each major subroutine are contained in the following sections. To 
retain flexibility and growth potential, the program modules are designed 
with minimum interdependence. 

Most of the technical and software experience used in designing 
FRACAS has heen obtained from work with STEAP (Reference 1) and the low 
thrust programs T0PSEP (Targeting /Op timization) , G0DSEP (linear error 
analysis), and SIMSEP (trajectory simulation). Many other programs were 
used in a lesser, but still significant, degree; POST (Shuttle trajectory 



optimization) , SWEAT (Swingby error analysis) , and BANANA (Bit Allocation 
Necessary for Accurate Navigation Analysis) . All of these programs were 
developed by Martin Marietta Aerospace and have been applied in a variety 
of interplanetary mission analyses. Some of the major technical analyses 
which were performed to develop algorithms are summarized in the Appendices. 

The appendices are self-contained memoranda complete with their own 
references. The two most difficult technical problems were in determining 
1) numerical accuracy of the covariance formulation and 2) method of covariance 
propagation including process noise model. These two problems were resolved 
satisfactorily and study results are summarized in Appendices 9.2 and 9.3 ? 
respectively . 


3 



2 . NOMENCLATUURE 


The following symbols are used throughout the program and subroutine 
descriptions. However, deviations from these symbols may occur in 
localized discussion if required for purposes of clarity. 


SYMBOL 


a 

c 

C 

E 

F 

G 

H 


K 

m 

P 

P 

« 

Q 

Q 

r 

s 

S 

t 

T 


u 

U 

U , V , w 
o o o 


DEFINITION 

-r ‘i ■■ ■ 1 ■ ■ " - 

propulsive acceleration 
propulsive exhaust velocity 
cross covariance 
target error index 
net cost function 
performance gradient 

observation sensitivity matrix WRT state 

parameter 

filter gain matrix 

spacecraft mass 

covariance 

propulsive power a t 1 AU 
dynamic noise matrix 
thrust noise matrix 
spacecraft position 
solve-for parameters 

target sensitivity matrix WRT control parameters 
arbitrary time 

event time, - target variables, or thrust 
dynamic consider parameters 
control parameters 

a priori covariances on dynamic consider 
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measurement consider and ignore parameters . 
respectively 



SYMBOL 


v 

w 

X 

r 

n 


y 

o 


T 


<P 


^k+l’V 


X 

SUBSCRIPT 


(. ) 
( ) 
( ) 
( ) 


A 

B 

C 

k+l,k 


( ) 
( ) 
( ) 
( ) 
(. ) 


S 

v 

W 

x 



DEFINITION 

spacecraft velocity or measurement parameters 
ignore parameters 
spacecraft state 
guidance matrix 

propulsive efficiency or time-varying thrust 
error 

transition matrix of dynamic parameters 

gravitational constant 

standard deviation 

correlation time of thrust error 

transition matrix of augmented state 

state transition matrix from time t, to t, _ 

k k+1 

target variation matrix 
DEFINITION 
assumed covariance 
true covariance 
state control covariance 
matrix evaluated over time interval 

‘k to \ + i 

evaluated at time t , t , t , respectively 

O K. X 

solve-for parameter 
measurement consider parameters 
ignore parameters 
spacecraft state parameters 
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SUBSCRIPTS 


DEFINTION 


( ) 


xu 


MISCELLANEOUS 
OD 
WRT 
E[ ] 

( ) + 

c r 

STM 


cross terms of state and dynamic consider 
parameters 

DEFINTION 

orbit determination 
with respect to 
expected value operation 
post^event value 
pre-event value 
State transition matrix 
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3 . PROGRAM DESCRIPTION (FUNCTIONAL INPUT/ OUTPUT) 

FRACAS is a pre-mission design tool used for parametric studies of 
trajectory dispersions and their relationship with anticipated error sources. 
It is an itermediate step between early mission opportunity definition and 
precision real-time flight software. As such, the FRACAS design reflects 
a trade-off between computational speed vs. high numerical and modeling 
accuracy. The results of FRACAS are intended to provide 1) trajectory 
sensitivities to dynamic processes, 2) state estimation accuracies based 
upon an orbit determination (OD) algorithm and expected errors in the 
environment, spacecraft performance, and navigation system, 3) trajectory 
correction requirements in the form of AV and/or thrust control adjustments 
to return the trajectory to desired terminal conditions, and 4) probabilistic 
trajectory dispersions as a result of all significant dynamic, guidance and 
naivgation processes. 

FRACAS is divided into three modes which, represent a logical sequence 
of analysis. The first mode generates a reference trajectory consistent 
with dynamic constraints. The user defines the mission in terms of launch 
and target planet, propulsion mode, flight time, etc. and then provides an 
estimate of desired control variables in the form of initial conditions and 
thrust parameters. The cpntrpl parameters are varied within constraints 
such that the final trajectory meets all desired end conditions and maximizes 
spacecraft mass at the target. In addition to providing a reference 
mission for use in the next two program modes, information is available 
relating to trajectory sensitivities and non-linearities with respect to 
dynamic parameters. 
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The second FRACAS mode is a linear covariance analysis of the reference 
mission. Distributions of errors in the form of variances and covariances 
are applied in a probabilistic sense to the trajectory as a sequence of 
mission events is processed. The mission events include thrust switching, 
guidance correction and navigation measurement processing (OD) . Usually, 
the time history of two types of trajectory errors are of interest: 
deviations of the actual trajectory from the reference, called control 
error, and deviations of the estimated trajectory from the actual, called 
knowledge error. Guidance events also provide probabilistic uncertainties 
in control corrections required to remove trajectory error at the event 
time . 

A key assumption in the error analysis mode is linearity, that is, 
deviations about the reference trajectory behave in a linear fashion. 

The third FRACAS mode verifies this assumption, or at least defines 
regions of linearity. Discrete errors ( randomly sampled from input 
statistical distributions) are applied to a deterministic trajectory. 
Guidance maneuvers are explicitly performed. By repeating the mission 
simulation with varying error samples, a Monte Carlo analysis can be 
constructed which takes into account the significant trajectory non- 
linearities. The simulation mode is of course the most lengthy in 
computer time and should be used primarily to support the error analysis 
mode . 

Together, all three FRACAS modes provide the analyst with trajectory 
data necessary for proper spacecraft subsystem and mission design. The 
program is designed to be structurally simple and easy to use, yet maintain 
flexibility with respect to more sophisticated analysis by applying existing 
options or by program change. 
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Each mode will have its own namelist input although many of the variable 
names will be common to more than one namelist, e.g,, basic spacecraft 
parameters. Printout options, punched output, and tape read/write will 
be controlled by namelist variables. At the beginning of every FRACAS 
data deck will be either an alphanumeric label or an integer which will 
determine the program mode. 

3 .1 Targeting and Optimization Mode 

The targeting and optimization mode (TOM) generates a reference 
trajectory which is supplied as basic input to the error analysis and 
simulation modes. The primary purpose of TOM is to incorporate in this 
trajectory all of the desired flight characteristics for a particular 
interplanetary or near-Earth mission while optimizing the final spacecraft 
mass. Injection conditions, a thrusting time history, and other control 
parameters are found which accomplish this optimization and yet lead to 
the required target conditions. The target constraints may be the final 
spacecraft state (cartesian or B-plane coordinates), final orbital elements, 
radius of closest approach, or other mission specifications which are listed 
in the input entries later in this section. 

Trajectories for the targeting and optimization mode are propagated 
using an Encke method with a two-step, fourth order Nystrom numerical 
integrator. Basically, the Encke method has been chosen to avoid integrating 
the entire vehicle acceleration vector to high precision. Since the 
accelerations due to a low thrust engine are considered small compared 
to the gravitational accelerations of the primary body, perturbation 
techniques can be applied. Integration is confined to evaluation of the 
relatively small deviations from the reference Encke conic resulting in 
rapid trajectory propagation. Conic propagation follows methods outlined 
in Battin (Reference 2) . 

9 




The trajectory generation mode features a discrete parameter iteration 
algorithm which accommodates the nonlinear aspects of the low thrust problem. 

The algorithm is a modification of the POST Shuttle projected gradient 
method (Reference 3) and uses finite differencing techniques which compute 
performance and target sensitivities to control variations. These sensitivities 
direct the control selection to maximize the performance index while 
minimizing the target error index. The performance index is simply the 
value of the final spacecraft mass while the error index is the weighted 
sum of the squares of the target constraint errors. 

The manipulation of trajectories to satisfy mission requirements is 
managed in the targeting and optimization submodes. TOM consists of 
four submodes which represent successive stages of trajectory development. 

These submodes are : 

1. grid generation 

2. trajectory targeting 

3. a combination of trajectory targeting and optimization 

4. trajectory optmization 

Generally, these submodes are employed in order as listed above. 

However, any submode may be skipped or used individually if the proper 
control profile is available. For example., the. linearity of controls 
characteristic of near-Earth missions, permits immediate entry into the 
targeting submode., although the control profile may not be. extremely 
accurate. This is not the case of most interplanetary missions where 
nonlinearities require accurate control estimates to be input and the 
grid submode to be implemented. In situations where either an interplane- 
tary or near-Earth mission is nearly targeted the third submode may be 
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employed initially. This submode provides control corrections for 
optimization of the trajectory in addition to completing the targeting. 
Finally, any trajectory which meets the targeting constraints can be 
optimized directly without entering the other submodes. 

The grid generation submode is available to produce a number of trajec- 
tories which do not necessarily satisfy mission requirements but provide 
a range of trajectory solutions. Thus, the main purpose of the grid 
submode is to locate desireable control regions for further examination. 

In turn, each control is incremented a fixed amount while the remaining 
controls maintain their nominal values. A single low thrust trajectory 
is generated for each control change and the associated target error index 
is calculated. Then pairs of controls are incremented and the target 
error indices are computed from the resulting trajectories. Subsequently, 
contours of constant target error may be plotted in the control space so that 
some control regions can be eliminated from further consideration. Upon 
completion of the grid the trajectory generation mode is terminated and the 
program user must choose the best control profile to initialze targeting 
and optimization or to employ another grid approach. 

When the targeting and optimization submodes are entered, a nominal 
trajectory is propagated directly from the input parameters. A series 
of tests is performed to determine which submode-targeting, optimization or 
b oth-is to be executed. If the target error index is large, the submode 
will be exclusively targeting. However, a target error index smaller than 
some arbitrary value (set in input) will result in simultaneous targeting 
and optimization. Whenever the index is below a specified lower bound, 
the optimization algorithm will be executed. 

After the submode decision the basic projected gradient method is 
applied to the controls. The targeting sensitivity matrix S and performance 
gradient G, are first computed. Elements of the S matrix represent the 
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sensitivities of individual target parameters to changes in controls and 
are used for both targeting and optimization. Similarly, the elements of the 
G vector represent the sensitivity of the performance index to changes in 
controls although these elements are used only for optimization. A 
weighting matrix which amplifies or diminishes the effects of the chosen 
controls is then calculated. Applying the projected gradient algorithm 
(Section 5.2), a new control vector direction is established. The magnitude 
of the control vector is determined by computing trial trajectories which 
adopt control profiles that lie in the new control vector direction. The 
new control profile is simply a scalar multiple of this control vector 
such that the targeting error index is minimized and/or the performance 
index is maximized. If the optimization is complete (the values of the 
performance index have converged to a maximum) TOM is terminated. 

Otherwise, the submode decision is made again and the cycle is repeated. 

The speed of convergence for various missions depends largely on good control 
estimates. One method of computing control profiles and sizing system require- 
ments for input into TOM is to apply the QUICKTOP program (Reference 4) 
developed by NASA/AMES to define low thrust interplanetary mission opportunities. 
QUICKTOP is approximate, self-starting, and computationally quick. The 
resulting values of the mission parameters can easily be adapted for 
refined targeting and optimization in the trajectory generation mode. 

Near-Earth missions, on the other hand, require less accurate control 
estimates. The input can usually be estimated by simple analytical calcula- 
tions . 

The targeting and optimization mode input is entered in the namelist 
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$DATA and is read in the main subroutine TOM, A second namelist $SMAG 


is used only when the targeting sensitivity matrix and performance gradient 
are to be input instead of calculated in the first iteration. 

$DATA 

Nominal spacecraft parameters 
o initial mass 

o base power and power supply constants (e.g. nuclear decay rate) 
o thruster efficiency 
Nominal thrust controls 
o thrust phase duration 
o pointing angles 
o thrust level 
o attitude mode 

Trajectory and integration parameters 
o numerical integration accuracy level 
o initial spacecraft position and velocity 
o initial epoch 
o trajectory termination epoch 
o launch and target planets 

o array of codes of intermediate gravitational bodies to be considered 
o trajectory stopping conditions 
oo sphere of influence 
oo radius of closest approach 
oo radius of designated final orbit 
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o Control parameter codes chosen from the following list 
of available parameters: 

initial spacecraft position and velocity, 
exhaust velocity, 

nominal thrust controls (e.g. thrust phase duration, pointing, 
thrust level, attitude modes) 

Targeting parameters 

o desired target parameter values (any combination of the 
following parameters) 

oo final spacecraft position and/or velocity 
oo hyperbolic approach velocity 
oo B-plane coordinates 

oo time of arrival at sphere of influence 
oo radius of closest approach to target body 
oo time of closest approach 
oo orbital elements 
o target tolerances 
Submode parameters 

oo grid generation 
oo targeting and/or optimization 
oo nominal trajectory only 
o maximum number of iterations 
o number of control parameters 
o number of target conditions to be satified 
o maximum change allowed in performance, in one iteration 
o limit of normalized targeting error below which "targeting only" 
is discontinued 
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o 


percentage of targeting error to be corrected in first iteration 
o estimated radius of linearity region in control space 
o maximum value of control change scale factor 
o curve fitting tolerances for trial trajectories 
o control parameter perturbations 
o control parameter weightings 

o minimum angle between control vector elements in the control 
space below which an element is deleted from control profile 

$SMAG 

o targeting sensitivities 
o performance gradient 

The trajectory information may be printed as a brief summary after 
each iteration or very detailed after each step of the projected gradient 
search. The detailed printout includes target sensitivities and weightings , 
performance gradients, trial trajectories, and control change scaling 
in addition to the desired spacecraft information throughout the optimized 
trajectory. 

3.2 Error Analysis Mode 

The error analysis mode performs a linear covariance analysis of 
guidance and navigation errors for low thrust trajectories. The under- 
lying assumption is that all trajectory errors may be described as 
linear deviations from a reference trajectory, and that their ensemble 
statistics are Gaussian. Verification of this assumption for any trajectory 
may be made by exercising FRACAS simulation mode, described in Section 3.3 
of this document. 

Probabalistic _a priori errors in the environment and spacecraft and 
tracking systems are propagated in time along the reference trajectory^ through 
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sequential events such as orbit determination (OD) and guidance corrections. 
Two types of ensemble error or covariances are distinguished - knowledge , 
which reflects the ability of the OD algorithm to estimate the spacecraft 
state; and control, which represents the dispersions of the actual 
spacecraft trajectory about the reference. Both knowledge and control 
covariances are stored internally in full covariance form rather than 
covariance square root form } the justification for which may be found in 
Appendix 9,2. Covariance propagation is done by either integration of 
covariance variational equations, or by the state transition matrix method. 

In general, the latter is recommended for reduced computer time (see 
Appendix 9.3). 

Error analysis flow proceeds sequentially from start time to each 
specified trajectory event. Event types availabe are measurement, 
propagation, eigenvector, prediction, thrust on/off, and guidance. A 
measurement event processes tracking data at a time point by applying the 
user specified OD algorithm. Available to the user are both Kalman- 
Schmidt (K-S) and sequential weighted least squares (WLS) filters. The 
filters are distinguished by their methods of gain matrix calculation. 

FRACAS modularity also allows the user to insert his own filter algorithm 
quite easily. 

A propagation event merely updates the knowledge covariance to the 
event time. Its primary value is in maintaining accurate covariance 
values during long propagations by forcing computation of the effective 
process noise over predetermined, user-specified intervals. Printout for the 
propagation event consists of the process noise covariance over the inverval, 
which may be suppressed at the user's option. 
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An eigenvector event converts all covariance matrix sub-blocks to 
variable standard deviations and correlation coefficients, all of which 
are output. It also computes eigenvalues, their square roots, and eigen- 
vectors for the position and velocity 3x3 sub-blocks of the state covariance 
matrix. Thrust on/off events are simply eigenvector events at the nominal 
thrust switching points. 

A guidance event is an update of the control covariance to reflect 
implementation of a trajectory correction. A correction is not performed 
deterministically, but only in a probablistic sense. Either impulsive AV 
or low thrust guidance can be performed (see Section 5.3.6). Low thrust 
guidance is further distinguished by being either primary or vernier. 

Vernier guidance is an update of a primary guidance correction to account 
for trajectory estimation improvement from tracking during the primary 
guidance interval. The guidance event computes, and displays to the user, 
expected correction covariances (Av or thrust control) , target error 
covariances before and after the guidance event, and the updated state 
control covariance. 

The simplest form of the filtering algorithm available estimates the 
six-dimensional spacecraft state - three position and three velocity 
components. Since there are always additional parameters whose uncertainties 
are important to the OD process, the error analysis mode is designed to 
accommodate these. Parameters may be included in two categories - solve-for 
parameters, which are estimated simultaneously- with the basic spacecraft 
state, and consider parameters, whose uncertainties are acknowledged by 
the filter, but which are not estimated. Consider parameters are divided 
into two types - measurement parameters which affect the measurement but 
not the dynamics, and dynamic parameters, those parameters which affect 
the dynamics and may or may not affect the measurements. 
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A major feature of the program is the inclusion of the generalized 
covariance option, a useful tool for studying filter sensitivity to 
mismodeling of real world error sources. When generalized covariance is 
exercised, two sets of knowledge covariances are operated on by the program. 
The first set, called assumed knowledge, comprises those covariances 
generated by the user selected filtering algorithm. The second set, called 
true knowledge, represents the effect the filtering algorithm has on true 
state estimation when real world error sources are not the same as those 
assumed by the filter. Mismatches between the two are effected either 
by setting true a. p riori uncertainties at different levels from assumed 
values, or the true state may be augmented by a vector of ignore parameters — 
parameters whose uncertainties are recognized by the true covariance analysis , 
but which are ignored by the assumed filter analysis. True covariance 
propagation and measurement updating is explained in Section 5.3.4, 
where the subroutine Filter is described. 

The most significant time saving option available to the user is the 
creation of a state transition matrix (STM) file. Since many different 
studies are often made on the same reference trajectory, the user may 
specify an event schedule which will include all time points at which 
events may occur. The trajectory generation overlay will then generate the 
state transition matrices between these event times and store them on 
tape. During execution of the error analysis mode, this STM file is read 
to retrieve the necessary transition matrices. If event times exist on 
the STM file between any two events in a specific error analysis, the 
transition matrices are multiplied together to compute the total transition 
matrix over that time interval. The use of the STM file considerably 
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reduces integration time for multiple studies of a single reference tra- 
jectory. For maximum efficiency in this multiple study usage, the genera- 
tion of the STM file must include transition matrix entries for all parameters 
which the user may at some time wish to solve^-for, consider, or ignore. 

When the error analysis recovers these matrices from the STM file, 
entries corresponding to current parameters are loaded into the proper 
transition matrix partitions, and those for unused parameters are passed 
over. 

Input to the FRACAS error analysis work is by namelists and event 
schedules where necessary. The first namelist ERRCON includes flags to 
indicate 1) if an STM file is to be created; 2) if the current run is 
supposed to execute an error analysis; and 3) if an error analysis is 
executed, whether the covaraince propagation is to be by transition matrices 
or integration of covariance variational equations. If either an STM 
file is to be created, or the covariance variational equations option 
is selected, the namelist ERTRAJ is required which includes all input 
needed for reference trajectory generation by either method. Since 
trajectory integration is required for prediction and guidance events, 
even when an existing STM file is used, all of the information in namelist 
ERTRAJ is written at the beginning of the STM file and is read from that 
file rather than cards. This guarantees consistency of integration 
accuracy level, gravitating bodies used, and nominal spacecraft control 
policy between the STM file and these event integrations. 

Immediately following ERTRAJ is a set of event scheduling cards defining 
all time points which must be written on the STM file. These cards are 
unnecessary if ERTRAJ is being read to initialize integration of covariance 
variational equations. 


19 



Next comes namelist ERANAL - which contains the basic information necessary 
for error analysis - followed by schedule cards for measurements and pro- 
pagation events. Last, if generalized covariance is to be used, comes 
namelist GENCOV, which initializes relevant parameters. Inputs to 
GENCOV are minimized by assuming that all true covariance information is 
the same as that for the assumed filter analysis unless changed by namelist 
GENCOV. 

Following are the error analysis mode namelists, and the input available 

through each: 

Namelist $ERRCON 

o STM file creation - true or false 
o Error analysis execution r- true of false 

o Integration of covariance variational equations - true or false 
Namelist $ERTRAJ 

o Initial spacecraft state and flag indicating coordinate system 
o Spacecraft mass, exhaust velocity, thruster efficiency 
o Flag indicating power source 

o Base power and power system constants, e.g. decay rate of nuclear power 
o Initial date 

o Final date or total flight time 

o Gravitating bodies to be used for trajectory generation 
o Integration accuracy level 
o Target body 

o Ephemeris of target body if not available internally 
o Parameter list for state transition matrices 

o Control array defining thrust on/off times, and nominal control 
policies for thrusting arcs 
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Namelist $ERANAL 

o All knowledge covariances describing augmented state 
o Parameter lists solve->for, consider, ignore 

o Time varying thrust parameters, their uncertainties and their 
correlation times 

o Filtering algorithm flag (K-S or WLS) 
o Control covariances 
o Print flags 

oo Print measurements according to time 
oo Print measurements according to type 
o Print measurements according to number, e.g. 

every 12th measurement 
oo Type of propagation event print 
o Number of measurement schedule cards to follow 
o Measurement noise levels 

o Station locations if additional or different stations from 
standard ones are desired 
o Event information 

oo Number of propagation event cards to follow namelist 
oo Number each of eigenvector, prediction and guidance events 
oo Event timing information 

oo Guidance policies and control weighting factors for each 
maneuver 
o Punch flags 

oo Knowledge and/or control punched at specified times to initialize 
later error analyses or for input to simulation mode 
oo Guidance variation matrices to eliminate recomputation in future 


error analyses 
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o Generalized covariance flag - true, or false 
Namelist $GENCOV 

o Ignore parameter list 

o A priori ignore parameter covariance terms 

o A p riori true covariance terms which differ from corresponding 
assumed terms 

o True time varying thrust parameter information which, differs 
from assumed 

o True measurement noise levels if different from assumed 
3.3 Simulation Mod e 

The purpose of the simulation mode is to examine trajectory non- 
linearities as they affect final target errors. Discrete a priori errors 
in the environment and spacecraft systems are applied as the trajectory 
simulation proceeds through each scheduled guidance and navigation event. 
The form of the simulation mode is such that many missions can be simulated 
quickly, each with varying samples of error sources, from which a Monte 
Carlo error analysis can be constructed. Tracking is simulated by sampling 
an estimation error covariance prior to each guidance event. Estimation 
error or knowledge covariances would be obtained from the results of 
linear error analysis. The sampled state error is added to the current 
actual state to form a best estimate which is used to design the maneuver. 
There are many options which can be selected for maneuver design, namely, 
choice of target variables , conditions and tolerances , linear or nonlinear 
(iterative) guidance, impulsive or low thrust corrections, thrust control 
parameter weighting and constraints, etc. After the maneuver is designed, 
execution takes place by applying the design maneuver plus execution 
errors to the actual trajectory. When all guidance events have been 
completed the actual trajectory is propagated to the target and end 
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conditions are evaluated. Statistical error characteristics for desired 
parameters are constructed after each, trajectory simulation. 

The sampling of estimation error covariances, as opposed to explicit 
orbit determination of the actual trajectory, was chosen because it was 
computationally faster, enabling a Monte Carlo error analysis to be a 
reasonable undertaking. We felt that Monte Carlo analysis provided 
much more information than a single trajectory simulation. The Monte 
Carlo approach also provides the flexibility of taking an "interesting" 
trajectory from the set of simulated missions and using it as a reference 
trajectory for linear error analysis. 

Because of the long run time necessary for a statistically significant 
Monte Carlo analysis, it is wise to break up the simultations into batches 
and keep the number of mission cycles per run to a minimum. Thus, capability 
exists in each FRACAS /simulation mode run to use the constructed error 
statistics of a previous run as- a- priori input and to punch cards containing 
cumulative statistics after the current run. 

Simulation mode input is divided into two namelists. The first 
($INSIM) is for describing the reference mission and associated errors. 

The second namelist (.$ INMAN) contains parameters describing a guidance 
correction event. Thus, each maneuver must have its own $INMAN. 

Namelist $INSIM 

o nominal spacecraft parameters; exhaust velocity, available 
thruster power, thruster efficiency, initial mass 
o variances in spacecraft parameters 
o nominal initial spacecraft state 
o state error covariance 
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o initial epoch 

o nominal thrust controls? phase duration, pointing., thrust level, 
attitude mode 

o thrust control variances (Bias) 

o time-varying thrust errors: mean and variance of correlation 
time, variance in thrust direction and proportionality 
o execution error variances for impulsive maneuvers 
o covariance of planetary ephemeris errors (.including gravitational 
constants) 

o launch planet, target planet, all other bodies to be considered 
o numerical integration accuracy level 
o random number initializer 
o print and punch flags 
o maximum number of mission cycles 
o number of maneuvers 

o number of mission cycles used to generate ji priori e rror 
statistics 

o cumulative a priori error statistics (from previous runs) 

Namelist $INMAN: 

o maneuver epoch 
o estimation error covariance 

o guidance law: linear or non-linear > AV or low thrust 
o guidance policy: cutoff condition (time, sphere-of-inf luence , 
closest approach, radius) and target set (B-plane coordinates, 
cartesian, conic, Earth synchronous) 
o target body 
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o target tolerances 

o thrust control tolerances and constraints 

o number of cycles used to generate- a pr iori error statistics 
o cumulative a p rior i error statistics (from previous runs) 
o target sensitivity or guidance matrix, target conditions, 
nominal spacecraft state and mass at maneuver epoch are all 
optional input 

Printout from the simulation mode can be a brief summary after 
each mission cycle or very detailed after each maneuver of each cycle. 
Cumulative statistics are always printed out at the end of the run and 
punched cards are available if desired. Some of the quantities displayed 
in the detailed printout will be deviations of actual parameters from 
their nominal values, cumulative means and standard deviations, target 
sensitivities and control steps after each iteration (for non-linear 
guidance), and reconstructed control and knowledge covariances. 



4. MACROLOGIC 


4.1 Functional F low 

The FRACAS program is a modular pre-flight analysis tool capable of 
generating targeted reference trajectories and performing error analyses 
on these trajectories as well as Monte Carlo trajectory simulations. 

FRACAS consists of three independent modules (T0M, TEAM and TSIM) 
illustrated in Figure 1. Each module performs the processing for its 
respective mode. FRACAS and the three modules are organized into an 
overlay structure to meet LRC imposed constraint of 70000 computer words on 

O 

CDC 6000 series computers. The program FRACAS is a main overlay while T0M, 
TEAM and TSIM are primary overlays. Because of extensive computational 
functions, TEAM is the only module which requires secondary overlays- . four 
in particular. Estimated core requirements for the entire overlay structure 
are illustrated in Figure 2. All estimates are based upon experience with 
MMA low thrust programs which perform functions similar to the proposed 
program. Should the estimates in Figure 2 be too optimistic . new 
secondary overlays could be created from the primary with subsequent increase 
in computer run time due to overlay loading . The total core requirement is 
estimated at about 67000 words. 

O 

A second LRC constraint of not more than 12 files has been met. 

FRACAS uses only four files; INPUT, 0UTPUT , PUNCH and a file (STM) used 
only in the error analysis module. 
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FRACAS 


Purpose: To generate a 

reference trajectory 
for low thrust missions 


Purpose: To perform a 

linear covariance 
analysis on a reference 
trajectory as a function 
of spacecraft, 
environmental, and 
navigational errors. 


Purpose: To perform Monte 

Carlo simulation of low 
thrust trajectories 
and related error 
s ounces . 


FIGURE 1. FRACAS Mode Structure 










Octal Thousands of Memory Locations 














4 .1.1 Targeting and Optimization Module (T0M ) 

T0M generates trajectories which satisfy specified target and control 
parameter constraints and maximizes the final spacecraft mass at the target 
planet. The module has the ability to operate in four independent 
sub-modes : 

o grid generation - finds desirable control regions 
o trajectory targeting and optimization - generates optimized and 
targeted trajectory 

o trajectory optimization ?- optimizes a targeted trajectory 
Generally, the submodes would be used in this sequence, however any 
subr-mode may be used if the proper control profile is available. 

The functional flow of T0M is illustrated in Figure 3. 

In the grid generation sub-mode a low- thrust trajectory is generated for 
each control change ('initial conditions and thrust parameters) and the 
associated target error index is computed. Then pairs of controls are 
changed and the target error indices are computed. Contours of constant 
error may be computed in the control space so that sections of the control 
region can be singled out for further study. 

In the targeting/optimization sub-modes a reference trajectory is 
generated or input which satisfies the control constraints and target 
conditions. If optimization is desired, changes are made to the controls 
to minimize a cost index. When a local minimum is found the optmization 
is completed. The projected gradient method is used for targeting and/or 
optimization. The targeting and optimization module is a single primary 
overlay' of FRACAS . 

4,1.2 Error Analysis Module (TEAM) 

The error analysis module is used to examine trajectory- dispersions 
resulting from thrusting, ephemeris , gravitational and measurement errors. 
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Figure 3. T0M Macrologic 


















Errors are represented as covariances which are propagated by state transi- 
tion matrices. Error covariances are updated by either a Kalman-Schmidt 
or sequential weighted least squares filter. The module is broken down 
into a primary overlay and four secondary overlays . The primary overlay 
contains all logic necessary to control initialization and cycling of the 
error analysis mode. The secondary overlays are defined as follows: 

1. DATA is responsible for all user input and editing; DATA will 
also do any initialization necessary for the proper functioning 
of the program. 

2. PATH generates the state vector and mass of the spacecraft and 
the transition matrices • from the previous to the current time. 

3. MEAS processes measurements by computing measurement noise and 
observation matrices, and updating the state covariance using the 
recursive estimation algorithm. 

4. GUIDM performs guidance events. 

Logic flow is shown in Figure 4. DATA is called to read the user's 
input and check for inconsistencies and omissions. DATA also performs 
some initialization such as zeroing variables and setting up event scheduling. 
For mulitple runs using the same reference trajectory, the user can 
create a file (S1H). containing the integrated state and transition 
matrices at each event time. On successive runs, the information from the 
file can be used instead of integrating the same trajectory repeatedly 
The basic cycle consists of obtaining the time of the next event ,• propagating 
the covariances to that time and calling the appropriate overlay to process 
the event, completing the cycle. 

4.1,3 Trajectory Simulation Module (TSIM) 

TSIM is used to examine the regions of linearity for low thrust 
trajectories as they affect target dispersions and thrust guidance 
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Figure 4. TEAM Macrologic 
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Figure 4. TEAM Macrologi:c (Continued \ 
































requirements. Errors are simulated by sampling error source distributions. 

The actual trajectory is propagated to a guidance or maneuver time where 
guidance corrections are designed using the estimated spacecraft state. The 
design maneuver is applied to the actual trajectory after execution errors 
have been added, and the actual trajectory is propagated to the next maneuver. 
When all maneuvers have been completed, the actual state is propagated to 
the target body and actual target conditions are computed. A Monte Carlo 
analysis is built from repeated passes through this basic cycle and statistical 
information is computed and printed. 

The module consists of a single primary overlay and is called only 
once for the entire simulation. The functional flow is shown in Figure 5. 

4.2 Subroutine Hierarchy 

As mentioned previously, FRACAS consists of three independent primary 
overlays or modes. The subroutine hierarchy for T0M, TEAM, and TSIM are shown 
in Figures 6 ,7 and 8, respectively, Multiple calls to subroutines are not 
shown but may be found in the detailed subroutine descriptions (Section 5) . 
Figure 9 illustrates the trajectory propagation hierarchy which is used in 
all three modes. Brief descriptions of these subroutines are given below 
along with references to detailed logic flow to be found in later sections. 


SUBROUTINE 

PURPOSE 

DETAILED 

DESCRIPTION (SECTION) 

BP LANE 

compute B-plane parameters 


5. 1.2. 2 

BUCKET 

sorts elements of a vector 


5.2.1 

C0yp 

controls propagation of covariances 


5.3.1 

CSAMP 

determines matrix eigenvectors/values and/or 
samples covariance 

5,4.1 

DATA 

processes error analysis input data 


5.3.2 

DATAS 

processes simulation input data 


5.4.2 
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SUBROUTINE 

PURPOSE DETAILED 

DESCRIPTION (SECTION) 

DELU 

computes control change 

5.2.2 

DERY 

computes covariance derivatives 

5. 1.2. 3 

DETECT 

detect changes in control 

5. 1.2. 4 

DYNO 

compute dynamic noise covariance matrix 

5.3.3 

ENCON 

compute or rectify reference conic 

5. 1.2. 5 

EP 

compute thrust parameters 

5. 1.2. 6 

EPHEM 

computes inertial state of a natural body 

5. 1.2. 7 

PEGS 

computes performance and target error indices 
and sensitivities 

5.2.3 

FILTER 

updates knowledge covariance by filtering 
equations 

5.3.4 

FUNCT 

selects trail steps for trajectory generation 

5.2.4 

GENMIN 

controls curve fitting for scale factor 
computation 

5.2.5 

GPRINT 

prints true estimation error statistics 

5.3.5 

GRAVFO 

computes gravity gradients and acceleration 
for primary and perturbing bodies 

5. 1.2. 8 

GRID 

generates grid of target errors in control 
space 

5.2.6 

GUIDM 

performs guidance events 

5.3,6 

GUIDS 

designs trajectory correction maneuver 

5.4.3 

INTEG 

performs 4th order Nystrom integration of 
R, V, and $ 

5. 1.2. 9 

MEN0 

computes measurement noise covariance matrix 

5.3.9 

N0ISE 

computes acceleration due to time-varying 
noise 

5,4.4 
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SUBROUTINE 


PURPOSE 


DETAILED 

DESCRIPTION (SECTION) 


PATH 

PINV 

P0WER 

PRED 

PRINT 

PR0P 

PTRAN 

RNUM 

SCHED 

SC0MP 

SETEVN 

SETUP 

SIZE 

STAPRL 

STAT 

STMGEN 

STMRDR 

TEAM 


controls reference trajectory generation: 
and state transition matrix computation as 
needed 

computes pseudo inverse of a matrix 
computes power output from low-thrust engine 
performs prediction events 
prints estimated error statistics 
propagates covariances 
computes STM 

generates a Gaussian random number 
determines time and type of next sequential 
event 

computes sensitivity matrix of target WRT 
thrust controls 

performs computation common to most events 
stores real-world or assumed constants into 
working arrays 

calculates magnitude of control change 
computes station location position and 
velocity partials 

computes cumulative statistics (mean and 

covariance) 

generates STM file 

reads STM file 

controls executive logic flow for error 


5.3.10 


5.1.2.10 

5.3.11 

5.3.12 

5.3.13 

5.3.14 

5.4.5 

5.3.15 

5.4.6 

5.3.16 

5.4.7 

5.2.8 

5.3.17 

5.4.8 

5.3.18 

5.3.19 
5.3.21 


analysis mode 



SUBROUTINE 


PURPOSE 


TEST 

T0M 

TRAJ 

TRAKM 

TSIM 

UPHILL 

USRGAN 

WEIGHT 

WLSGAN 

XGUID 


DETAILED 

DESCRIPTION (SECTION) 


tests for convergence 5.2.11 

controls 1/0 and initiates targeting/ 5.2.12 

optimization mode 

controls Encke integration 5. 1.2.1 

computes observation matrices 5.3.22 

controls logic flow for simulation mode 5.4.10 


contols logic flow for targeting/optimization 5.2.13 
computes filter gain matrix with user supplied 5.3.23 
algorithm 

computes weighting matrix 5.2.14 

computes filter gain matrix according to 5.3.24 

sequential weighted least squares algorithm 
controls execution sequence for guidance 5.3.25 

events 
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FIGURE 6 < 


Subroutine Hierarchy for Targeting and Optimization Mode 
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FIGURE 7. Subroutine Hierarchy for Error Analysis Mode 






































FIGURE 8. Subroutine hierarchy for Trajectory 
Simulation Mode 






































5 SUBROUTINE DESCRIPTIONS 
5,1 Subroutines Used in More than One Mode 


5,1.1 Utility Routines 

Many small routines are used by several of the FRACAS modules . 
These are described briefly below: 


o ADD 


t- add two matrices 


o ADPR - additive matrix product, [a] <= [a] + [Bj |cj 
o ADPRT - additive matrix product [X] = |Aj + [Bj [cj 

o ADXYXT - additive matrix product ;A] = | Aj + |b] [cj [B^ 

o C0PY - copy one matrix into another 

o C0PYT - copy the transpose of one matrix into another 


o C0RREL 


compute standard deviations and correlations of 


a covariance 


o DD0TB - compute the inner product of two vectors 
o DUMAG compute the magnitude of a vector 

o DXB - compute the cross product of two vectors 

o GHA compute the Greenwich hour angle 


o INVERT - invert a matrix 

o JAC0BI - compute the eigenvectors (eigenvalues of a matrix) 
o MAT0UT - print a matrix 

o MULT - matrix product [a] = I B] [c] 

o MULTT - matrix product [Aj = [BJ [cj T 

o PINV - pseudo inverse of a matrix 


o SDC0V 


converts standard deviations and correlations 


to covariance 


o SUB - subtract one matrix from another [Aj = [b] - [cj 

o SUBT - subtract one matrix from another. [A] = [b| - [C! 

o SYMTRZ - symmetrize a matrix 
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o TADPR 
o TADPRT 
0 TIM 

o TMULT 
o TMULT T 
o UNITV 
o XTYX 
o XYXT 
o ZMAT 


- • T r- 

- additive matrix product [a] = jAj + ;B] ,c] 

additive matrix product [a] = [a] + . B] [c] 

- convert time in seconds to days , hours , minutes , 
and seconds 

■"X r 

- matrix product jAj = [B] |c] 

r ">X X 

- matrix product [a] = [Bj [(^ 

- unitize a vector 

x 

- matrix triple product [Aj = [Bj [c :Bi 

. , r {E 

- matrix triple product ^ A1 = , Bj jCj [JBj 

^ zero out a matrix 
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5. 1.1.1 Subroutine; 

Purpose : 
Input : 

Output : 


Remarks : 


then 


C0NIC 

To convert cartesian coordinates to conic elements 
o position, r_ 
o velocity, v 

o gravitational constant of primary body, y 
o semi-major axis, a 
o eccentricity, e 
o inclination, i 
o longitude of ascending node.fi 
o argument of periapsis ,w 
o mean anomaly , M 
let h = r_ x v 
£ 

w = a 

d = j: • v 

c = — (vxh) - r 
- y 1 ' - 

P = c 

s = |h|/y 

i = cos "*"(w ) 

— z 

i w 

n -1 -X 

fi = tan — 

-w 

-y 

q = w x 

“ = tan_1 ( ^/q z } 
sin(9) = (|h|d)/|rj 

cos (9) = (|hj 2 * * * * * 8 -y ) / | rj 

9 = tan ^ (sin (9) /cos (9)) 

8 = cos 1 (d/(.|r| |v| ) 
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cos(e) = 1 

a 

sin(E) = d/ p |a| 

for elliptical case (a>0) 

r = tan ^ (sin (E) /cos (E)) 

M = E - sin(E) 
for hyperbolic case (a<0) 

sinh(f) = sin(E)/|cJ 
cosh(f) = cos(E)/|c_| 

E = ln( s inh(f) + cosh(f)) 
M = sin(E) - E 
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5. 1.1. 2 Subroutine 
Purpose : 

Input ; 

Output : 


C0NVRT 

To convert spherical coordinates to cartesian 
coordinates . 

o spherical coordinates of position (r , 0>9) 
o spherical coordinates of velocity (y y,o) 
cartesian position vector R 
cartesian position vector V 


Remarks : 


R 

x 

R 

y 

R 

z 

B 

x 

B 

y 

B 

z 

V 

x 

V 

y 

v 

z 


r cos0 cosG 
r cos0 sinG 
r sin0 
v siny 
v cosy sino 
v cosy coso 

B cos0 cosG - B sinG - B sin0 cosG 
x 'y 

B cos0 sinG + B cosG - B sin0 sinG 
x y z 

B sin0 + B_ cos0 
x *• 
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5. 1.1. 3 Subroutine EULMX 


Purpose : 

Input : 

Output : 
Remarks : 

where 


To compute a rotational transformation matrix from 
the Euler angles, 
o Euler angles (a,S, Y) 
o axes of rotation, a , i = 1,3 
o transformation matrix [p] 

[p]= [h] [g] [f] 

[f] = f(a, a^) 

[G]= f(3, a 2 ) 

[K]= f(Y, a 3 ) 
f(i|>, a) is defined as 



' 1 

0 

0 

1 

for a = 1^ f (xp , a) = 

0 

cosip 

sin^ 


. 0 

-sin^ 

cosiJj 


cosi|j 

0 

-sinipl 

for a = 2 ; f (if), a) = 

0 

1 

0 



. sinij; 

0 

COSlJj J 


cosij; 

sinijj 

0 “ 


for a = 3 ; f (jp , a) = 

-sinijj 

cosij; 

0 



. 0 

0 

1 j 
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5. 1.1. 4 Subroutine 


PECEQ 


Purpose : 


Input : 


Output : 
Remarks : 
Let 

A 

P = 

and 

A 

N = 


then 

where 


To compute the transformation matrix from 
planetoecentric ecliptic to planetocentric 
equatorial coordinates, 
o planet number 
o Julian date 

o planets conic elements, (a,e,i,n,w, M) 
o planets right ascension a and declination 6 of 
the pole vector 

o obliquity of the ecliptic, e 
transformation matrix [a] 


A 

P be the planetary pole vector, 
cos a cos 6 

cos e sin a cos 6 + sin e sin 6 

-sin e sin a cos 6 + cos e sin 6 

sin i sin R 

-sin i cos Q - 

cos i 

Tl r A < A 1 / ^ T 1 

[a] = [x ; Y j Z] 

A A 

Z = P 

X = P x N/|p x n| 

f t. f 


Y = Z x X 



5 . 1.2 
5 . 1 . 2.1 


Trajectory Routines 


Subroutine 
Purpose : 

Input : 


Output : 


TRAJ 

To control the integration of the trajectory (and certain 
other parameters) between two time points 
o thrust controls 
o true state _r,v 
o covariance integration flag 
o orimary body 
o target planet 
o start time, t^ 
o stop time, tk+i 

o dimension of state transition matrix, n 
o true state at t^ + ^ 

o integrated state transition matrix or augmented state 
covariance matrix 


Remarks : 

TRAJ is the logic control routine for the integrator. The Encke 
perturbed conic method (Ref. Battin ch. 6) is used with a Nystrom fourth- 
order two-step numerical integration technique. TRAJ can optionally 
integrate the covariances directly or compute state transition matrices 
for either the basic state or the state augmented by thrust parameters. 
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TRAJ-3 
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TRAJ-4 













































TRAJ-7 



RETURN 




5. 1.2, 3 Subroutine DERY 


Purpose ; 


Input : 


Output : 
Remarks : 

Let P = E 


To compute the derivative of the augmented state covariance 
matrix . 

o augmented s tate covariance matrix [ Pj 
o thrus-t transformation matrix [T] 
o gravity gradient [g] 

o process noise correlation times, t , i=l,6 
o augmented state covariance derivative matrix [Pj 



where X = 


= deviations of position components from nominal 

v = deviations of velocity components from nominal 

r\ = deviations of thrust components due to noise 

ii = deviations of thrust components due to bias 

T 


then P = [Fj [Pj + [P] [F] 


where P is partitioned 

P = 


and F is partitioned 

F = 



C P 

pv v 

C C 

pn vn 

c c 

pu vu 

0 I 

G 0 

0 0 

0 0 


T 

C 

pn 



P 

n 

c 

nu 

0 

N 

H 

0 



0 

T 

0 

0 
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DERYt-2 


where 

N = 

[g g] 





H = 


0 





0 

- 1/X 6 



then 

Pp] 

■ i C pv] 

+ [ C pv] T 



[ C pv] 

- [ G U 

p 1 + 

pj 

r n] r c ] 

! J i ppJ 

+ T 1 " C ] + I p 1 
u J * pu J 1 V J 


[ C p„J 

‘ [=][ 

c ! + r c ’ 

pnj vn - 



IV] 

= [C 1 
L vuj 





PvP 

: g h c pvJ t 

+ l t j[ 

G J + [ C pJ[ G J T + Kn] T ( N ] T + 


[M 

- f H "i [" C | 
L J L v J 

+ t c , 

Ji G J T+ [ p 

ni-" N J T + f C nu] T i" G ] T 


re 1 

1 . vuj 

= 'C If G 
i puJL 

; T + 

[ c nu i[ N J 1 

+ [ P u)[ G J T 


FP 1= 0 
L nJ 

W M'nJ 

l>u] ' 0 
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5. 1.2. 4 Subroutine 


DETECT 


Purpose To detect control changes- during an integration step 

and break up the step at the time of change 

Input: • current time t 

• proposed step size h 

• control times T 

c 

Output : None 

Logic flow: 










5. 1.2. 5 Subroutine ENC0N 


Purposet To propagate reference conic to current time or 

rectify conic if deviations are too large. 
Input: • true position vector 

• true velocity vector 

• osculating conic 

• previous time 

Output; • updated osculating position vector 

• updated osculating velocity vector 

• osculating 



61 








5. 1.2. 6 Subroutine 


EP 


Purpose; To compute magnitude and direction of low thrust. 

Input; o spacecraft mass, M 

o exhaust velocity, c 
o base power, P^ 
o engine efficiency, n 
o spacecraft position, R 
o spacecraft velocity, V 

o thrust controls: control type, thrust scale factor(s) 
o simulation flag 

o simulation error levels, 6E_^, i=l,3 
Output; o thrust vector, y 

o thrust vector rotation matrix into inertial 
coordinates, [l] 
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EE-2 
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Cone-clock 

control 


a = a * s 

Compute current cone-clock angles 


Discrete thrust errors 


a = a (1 + (jEi) 

cone = cone angle + <jEg 

clock = clock angle + <JE 3 


X = a 


COS (Clock) SIN (Cone) 
SIN (Clock) SIN (Cone) 
COS (Cone) 


Compute right ascension a and declination 5 of R in inertial 
coordinates , 


COSa SIN <5 
SINa SIN(5 
- cos <5 


sin a cosfl cos a cos 6 
cos a cos$ sin a cos 6 
o sin 6 


= [A] X 











a 


.5a 



Compute 

initial 

right ascension a and declination i 
coordinates 

of Y in 

[T] = 

cos a COS 6 
SIN a COS 5 

-SIN a cosi 
cos a cosi 

-COSa 
sin a 

SIN 6 1 
SINi 


SINa 

0 

cos 5 J 









5. 1.2. 7 Subroutine 


EPHEM 


Purpose : 
Input : 


Output ; 
Logic flow: 


To calculate the position and velocity of a planet 

• Julian date 

• planet code(s) 

• planetary constants 
heliocentric state vector of planet 
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5. 1.2. 8 Subroutine 
Purpose : 

Input ; 


Output ; 


ORAVF0 

To compute gravity gradient o£ primary body and 
perturbing bodies- and to compute accelerations 
caused by low thrust propulsion. 

o osculating spacecraft position vector relative 

to nrimary body- r 

— osc 

o difference between osculating position and 
true position, 6_ 

o planets or bodies to be considered, B(i), 
i=l,...,N 

o positions of perturbing bodies relative to 
primary, p. 

o gravity gradient matrix, G 

o acceleration due to perturbing bodies and 

thrust, a, 
d 
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Logic Flow: 


GRAVF0-2 
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5 . 1 . 2 . 9 Subroutine INTEG 


Purpose : 


Input : 


To numerically- integrate the equations of motion 
C and the variational equations if desired) over 
an integration step. 

o initial state deviations from conic r , v 

o ' o 

o perturbing accelerations, 

o state transition matrix at start, 

o gravity gradient matrix, [g] 

o current spacecraft mass, m 
o exhaust velocity^ c 
o propulsive efficiency, q 

o step size, h. 

•A — X 

o true state, R,V 
o covariance integration flag 

o augmented state covariance [p] 
o augmented state covariance derivative [p] 
o Thrust controls 
o thrust acceleration, T 

o dimension of state transition matrix, n 
o thrust transformation matrix [f] 
o current time , t 

o osculating spacecraft state r , v 

r -osc -osc 

2 

o mass variance , a 

M 

o acceleration proportionality variance, a^ 

o correlation time , x 

o acceleration scale factor, a 

* s 

o acceleration resolution variance, a 

r 
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Output ; 


o 


o 

o 

o 


integrated state deviations r r , v. 

— x — f 

updated true state, R,V 
2 

mass variance, a,, 

M 


updated state 


transition matrix 


partitions 



o updated augmented covariance matrix, [p 
o current time, t 


Remarks j 

The numerical integration technique is a fourth order Nystrom. 
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INTEG-7 
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INTEG-9 














5.1.2.10 Function 


P0WER 


Purpose : 


Input ; 


To compute the. power ratio available to the propulsion 
system. 

• model selection 

• flight time, t 

• heliocentric distance (_for solar propulsion) ,r 


• power constants, 

• range of usefulness for solar array, 


(P/P o ) 

max 


Output: • power ratio, P/P^ 

Logic flow: 
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RETURN 












5 .2 Targeting Optimization Mod e 


^•2.1 Subroutine 
Purpose : 


Input ; 


Output : 


BUCKET 

To sort a set of independent elements in ascending 
order and to find a bounded minimum from the 
associated set of dependent elements, 
o set of independent elements, X^. 
o set of dependent elements. Y_^ 
o number of elements, N 

o ordered set of independent elements, X. 
o ordered set of independent elements, Y. 


o pointer, k^ to a minimum dependent element 

Remarks : 

This routine is used in preparation for the polynomial curve fitting 
routine, MINMUM, to aid in calculating the new control profile. 

BUCKET sorts pairs of elements (X^, Y ) in ascending order of the 
elements X^ and locates the element Y^ from the newly ordered pairs such 


that 


Y <Y <Y 
k-1 k k+1 

If this condition cannot be satisfied the pointer, K, is set to zero to 
indicate that no bounded minimum exists. 
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BUCKET- 2 














5.2.2 Subroutine 
Purpose ; 
Input s 


Output : 


DELU 

To compute the control correction vector, 
o snhmode designation 
o sensitivity weighting matrix, W 
o target sensitivity matrix, S 
o performance gradient, G_ 
o target errors, A^T 
o current control vector, U 
o estimated radius of region of linearity 
o number of controls, M 
o number of targets , N 

o complete control correction vector AU_ 
o optimization control correction vector, AU^ 
o constraint control correction vector, AIL-, 


Remarks : 

Subroutine DELU applies the projected gradient algorithm to compute 
the control correction vector, All. The direction of the correction vector is 
dependent upon the submode designation. For example, 

Targeting only: A^U = AU^ 

Targeting and optimization: AU_ = AU^ + AIJ^ 

Optimization only: Al[ = AU^ 

Linearly dependent controls are identified in subroutine STEST and 
are dropped from the subsequent matrix operations. No change is allowed 
in the omitted controls for the current iteration. 
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Logic flow: 


ENTER 


DELU- 2 
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5.2.3 Subroutine 
Purpose : 


FEGS 

To calculate the performance index, the error index, 
the targeting sensitivity matrix, and the performance 


Input : 


Output : 


gradient , 

o desired target values, T 

' —o 

o number of targets, N 
o nominal controls, U 
o number of controls, M 
o control perturbations, 6U_ 
o flag to indicate desired computations 
oo generate nominal trajectory only 
oo compute performance gradient G_ and target 
sensitivity matrix S only 
oo generate nominal trajectory, (1 and S 
o performance index, F 
o target error index, E 

o values of target parameters for nominal trajectory 
o performance gradient , (1 
o targeting sensitivity matrix, S 


Remarks : 

The performance and target error indices which are computed in 
FEGS are used in subroutine TEST (section 5.2.11) to determine the 
routine submode for the next iteration. The performance index is simply the 
final spacecraft mass and the error index is the sum of the squares of the 
target error. 
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Logic flow: 


ENTER 


PEGS -2 












FEGS -3 



i 







5.2.4 Subroutine 
Purpose : 

Input : 


Output : 


FUNCT 

To calculate the net cost-function for a trial 
trajectory . 

o current control vector, U 
o trial control change scale factor, y 
o control change vector. AU 
o current performance index. F 
o current sensitivity matrix, S 
o current performance gradient, G 
o desired target values, 
o submode designation 
oo targeting only 
oo targeting and optimization 
oo optimization only 

o net cost-function value , F Cy) for the trial 
trajectory 


Remarks : 


The net cost-function is described in Section 5.2.8 (Subroutine SIZE) 
F (y) = aF T (y) + g F Q (y) 


a = 


1 for targeting only or simultaneous targeting and optimization 
0 for optimization only 


g =< 


jl for optimization or simultaneous targeting and optimization 
[o for targeting only 
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Logic Flow 


ENTER 


FUNCT-2 



1 

V 


9 . 








FUNCT-3 





5.2.5 Subroutine. 
Purpose : 


Input : 


Output ; 


GENMIN 

To generate a series of trial trajectories based on 

control change -vectors of different magnitude and 

to choose the best control change scale factor. 

o current net cost- function value, F(y) 

■ y=0 

o value of the first derivative of the net cost- 

t 

function evaluated at y-0 , F (0) 
o curve fitting tolerance for trail steps, n 
o maximum value of y 

o value of the net cost-function for each trial 
trajectory 

■k 

o minimum value of net cost function, F(y ) 

: k 

o minimizing scale f actor , y 


Remarks : 

The net cost-function is described in Section 5.2.8 (SIZE) 
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GENMIN-3 



? 








GENMIN 



Approximate F(y) with a second order Poly- 
nominal (P) ; coefficients based on 
three sample points. Compute y and 

P 3 (y 3 *> 








GENMIN-5 
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5.2.6 Subroutine GRID 

Purpose: To generate a family of trajectories. 

Input: o nominal controls, JJ 

o control increment, AIJ 
o number of controls, M 

o maximum value of scale factor y for control 

max 

increment 

o desired target values 

o flag designating two incremented controls per 
grid trajectory 
Output : o External 

Remarks : 

Subroutine GRID is used to generate a grid of trajectory target 
error indices. The error indices are used either to direct a finer grid 
search or to choose a control profile to enter the other submodes of TOM. 
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5.2.7 Subroutine 


MINMUM 


Purpose : 


Input : 


Output 5 


To find the minimum value of a function, F(y) : , and the 
minimizing independent variable, y 
o flag denoting type of polynomial approximation 

oo second order polynomial, coefficients based on 

two sample points and one derivative evaluated at 
a point 

oo third order polynomial, coefficients based on three 
sample points and one derivative evaluated at a point 
oo second order polynomial ^ coefficients based on three 
sample points 

o set of at most three distinct values of the function, F(y) 
o set of corresponding independent variable values , V 
o value of the first derivative of the function F (y) 
evaluated at y=0 

o estimate of the minimum value of the function ,F (y) 

"k 

o value of the minimizing parameter y 


Remarks : 

The function, F(y), is approximated by either a second or third 
order polynomial, P(y), in order to compute analytically the minimizing 
parameter y . The polynomial approximation is of the form 


F(y) = P (y) 



i=0 


V 


i 


where n— 2 or n-3. The following three cases describe the method of approxima- 
tion and the resulting minimization process. 
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MINMUM-2 


Case 1 


F is fitted with a quadratic polynomial based on; 

1) F (0) 


2) F (0) = 


dF (y ) 
dy 


y=0 


3) F(y ) where y >0 is an initial estimate of v 
o o ' 

The quadratic polynomial coefficients are calculated from 
the formulae 
a = F (0) 


a x = F (0) 


a„= F (y ) 7 

2 o Y 


a a 

0 +JL 


The independent variable value minimizing the quadratic is 


* 

y = 


_L 

2a. 


Case 2 F is fitted with a cubic polynomial based on:. 

1) F (0) 

2) F CO) 

3) F (y ) where y is as in Case 1 

o o 

4) F(y^) where y^>0 is a sample value 

Tfie cubic polynomial coefficients are calculated from the following formulae 
a Q = F(0) 

a 1 >= f'(0) 


a 9 >= F(qA) -a F(A ) -Aa(l+a)a 
(1-a) 1 


(l+ a+ a ^) S q 


Xa.a+ a (1+a) + 
a. - 1 o 

J m - 

A 3 cr 


F(a) - F(aX) 

■Z kk 
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MINMUM-3 


where X = max (y .y ) 

o ' 1 

= min (A ,y )/X 
o i 

* 

The independent variable value, y , minimizing P is 

y* = ( ~ a 2 + a 2 2 " 3a 3 a i ) 

3a„ 


Case 3 


A quadratic polynomial is fitted to F(y ) s F(y ) and F(y.) 

2 3 4 


where Y 2 ,Y 3 - anc * ^4 are § reater than zero and represent sample values of y, 
It is assumed that the input of satisfy two conditions 

1) y 2 <y 3 <y 4 

2) F(y 2 ) F(y 3 ) F(y 4 ) 

The formulae for the quadratic coefficients are as follows; 

b . . = y.y . 
ij i J 

c. . = y . + y . 
ij i J 

d. . = y . - y . 
ij i J 


a = 


° d 23 d 24 


b i4 b ?A b-y-i 

' F( V + 1— sr F( V + rf: F <V 


32 34 


42 43 


a l d„„ d 


F (y ) F (y ) F (y ) 

V d„ d„, d,„ d_ n V 


23 24 


32 34 


42 43 


f(y 2 ) F(y 3 ) f(y 4 ) 

a 2 d„„ d„ , d„„ d„ , + d,„ d, 


23 24 


32 34 42 43 


The independet variable value is the same as in Case 1. 


Y = - 


2a. 
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Logic flow: 


, ENTER 





5 












5.2.8 Subroutine 


SIZE 


Purpose: To calculate the magnitude of the control change 

vector. 

Input: o current control vector, IJ 

o value of the performance index, F 
o percentage of target error to be corrected in 
one iteration 

o values of target errors , AT 
o target tolerances 
o target sensitivities, S 
o performance gradient, (5 
o estimated size of region of linearity 
o submode designation 
o individual control scale factors 
o type of weighting matrix to be computed 
o initial estimate of control change scaling factor, y 
o curve fitting tolerance for trial steps, n 

Output : o complete control change vector , AU_ =AU^+AIJ 2 

o optimization test angle, 0 

Remarks : 

Prior to calling subroutine DELU, the targeting error correction 
for the current iteration is computed. The nonlinear effects of 
certain targeting problems require that only a certain percentage of the 
target error is removed in any one step to prevent divergence. 

For any particular control vector- IQ in the independent-variable 
(control) space the projected gradient algorithm reduces the multi- 
dimensional problem to a one dimensional search either along the constraint 
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SIZE-2 


direction to minimize the sum of the squares of the constraint violations 
or along the optimization direction to minimize the estimated net cost- 
function. In either case, once the initial control vector U_ and the 
direction of search AU are specified, the problem reduces to the numerical 
minimization of a function of a single variable - namely the scaling 
factor y • 

Subroutine GENMIN is called to compute the value of the scaling 
factor y * which minimizes a function F(y) in both the constraint 
direction (AU^) and optimization direction (AU^)or each direction 
individually depending on the submode designation. The net cost- 

function is the sum of two functions, F fy) and F (y) . 

T o 

f(y) = <*f t (y) + £ f q (y) 


OL - 


1 for targeting only or simultaneous targeting and optimization 
0 for optimization only 


3 = 


1 for optimization or simultaneous targeting and optimization 


0 for targeting only 


The first derivative F (y) 


Y=0 


is used in the one dimensional search to 


find y and is calculated in SIZE prior to the call to GENMIN. 

The function F^(2f) to be minimized along the constraint direction 
A- U 2 is the sum of the squares of the target errors 
f t (y) = | |atl(u+yau 2 ) 1 1 2 

The -first derivative evaluated at » =0 is then 


' T 

F (0 ) = 2AT (U)SAU, 
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SIZE-3 


The function F q ( TO to be minimized along the optimization directional]^ 

is the estimated net cost-function which is defined 

F o (y) = FOJ+yAUj) -F(U) + G T (U) [-S T (SS T )" 1 AT(U-b Y AU_ 1 )] 
v . ' ' 

Change in performance Linearized approximation to change 

index produced by a step in performance index required to 

of length y along AU^. maintain the current target errors. 

The first derivative evaluated at y= 0 is then 

f' (0) = g t (u) AU. 

o — — ”1 

Hence 

F* CO) = of’(0) + sf'(o) 

i o 
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S£ZE-4 
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SIZE-5 
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5.2.9 Subroutine 


STEP 


Purpose 
Input ; 


Output : 
Remarks 
The new 


To compute the new control vector 

o control vector. U , , 

—old 

o control vector scale factor, y 

o control vector change AU_ 

o dimension of the control vector, M 

o new control vector U 

—new 

control step is U = U , , +yAU 
—new —old — 



5.2.10 Subroutine STEST 


Purpose : 


Input ; 


Output ; 


Remarks ? 


To determine linearly dependent controls among the 
elements of the control correction vector 
o sensitivity matrix (the partial derivatives of the 
target variables with respect to the controls) 
o number of target variables, N 
o number of control variables, M 

o tolerance value, e, determining linear dependency 
between two controls 

o number of linearly dependent controls 
o those controls which have been eliminated from the 
control profile as linearly dependent 
The inner products between the columns of the sensitivity 
matrix , S , are computed where 


S = 


3T ] 

3u; 


3T 2 

3^ 


3T 

3tT 


3T ] 
’ 3U~ 


m 


3T 

3U 


N 

1 



S1 x 
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S TEST-2 


If the value of the inner product, P, between two columns, I and J, is 
such that 

(1-e) <_ P <_ 1 

then the controls, U^. and U , are considered linearly dependent and one of 
the controls is eliminated from the control profile for at least one 
iteration (there exists the possibility that within a different region of 
the control space, P will not satisfy the preceding test condition and 
the control may again be added to the control profile) . For any given 
pair of linearly dependent controls, the first control is arbitrarily 
eliminated from the profile unless the second control appears in one or 
more other linearly dependent pairs. If this situation occurs the second 
control is eliminated from the profile for at least one iteration. 
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Logi c Flow; 


STESTt-3 



RETURN 









5.2.11 Subroutine 
Purpose : 

Input ; 


Output : 


TEST 

To test for convergence and to make a decision for 
targeting and/or optimization in the next iteration, 
o iteration number 
o maximum number of iterations 
o target error index, E 

o limit to which E must be reduced before the 
targeting submode is discontinued, T 
o lower bound of E below which the simultaneous 

targeting and optimization submode is discontinued,!^^ 
o optimization convergence test angle, 9 
o angle below which the optimization is considered 
complete , e 
o submode flag 
oo target only 
oo target and optimize 
oo optimize only 
o convergence flag 

oo iteration converged 

oo iteration not converged 

oo maximum number of iterations reached 


Remarks ; 

The iteration is considered converged and the run is terminated 
when the performance index is maximized. The test angle 9, which approaches 
zero as the. optimization is completed, is a means of testing the convergence 
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TEST-2 


status (Section 5.2.8 - SIZE). 

The decision for targeting and/or optimization is based on the 

current value of the targeting error index, E. If the value of E is 

greater than T the targeting submode will be entered. If 

T, <E<T 
low up 

then simultaneous targeting and optimization will occur. A value of E 
less than T^ will result in optimization only. 
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5.2.12 Subroutine TOM 


Purpose: To initialize all parameters and to choose the 

proper trajectory generation algorithm. 

Input: o see the input description of Section 3.1 

(Functional Input/Output) 

Output: o see output description of Section 3.1 

(Functional Input/Output) 

Remarks : 

TOM initializes parameters through the namelist input $DATA. 

The necessary parameters which are not input assume the default values 
set in TOM. If the targeting sensitivity matrix and performance gradient 
for the first iteration are also to be input the namelist $SMAG will be 
read. 
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T0M-2 
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5.2.13. Subroutine UPHILL 


Purpose : 


Input ; 


Output : 


To generate a targeted and optimized reference 
trajectory. 

o number of target constraints 
o number of controls (independent variables) 
o maximum number of iterations 
o initial estimate of control vector 
o error tolerances for targets, € 
o upper RMS constraint error tolerance 
o lower RMS constraint error tolerance 
o estimated radius of region of linearity 
o percentage of error to be corrected on first 
iteration 

o sensitivity matrix S and performance gradient G 
for first iteration 
o iteration number 
o value of performance index 
o value of the error index , E 
o optimal control vector, U 
o control change scale factor, y 
o total control correction ,Al[ 
o optimization control correction ,AU^ 
o constraint control correction 
o target sensitivities, S 
o performance gradient , G 
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UPHILL -2 
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5.2.14 Subroutine WEIGHT 

Purpose: To generate a weighting matrix 

Input: o nominal control values, U 

o number of controls , M 
o targeting sensitivity matrix, S 
o desired target tolerances, £ 
o number of targets, N 
o flag designating type of weighting 
oo control weighting 
oo sensitivity weighting 
Output: o diagonal weighting matrix, W 

Remarks : 

The weighting matrix is used in the projected gradient algorithm 
to emphasize other controls. Two options are available: 1) W based on 

'i 

the largest modulus of the sensitivity elements of each row of the S 
matrix (see Section 5 . 2 . 10^STEST ; f or a description of S) , and 2) W 
based on the square of the control values. W is a diagonal matrix with 
all off-diagonal terms equal to zero. 
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WEIGHT -2 


Logic Flow 
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5.3 


Error Analysis Mode. 


5.3.1 


Subroutine 
Purpose : 

Input : 


Output : 


C0VP 


To propagate a covariance matrix from one event time 
to another 


o 

o 

o 

o 

o 

o 

o 

o 

o 


initial epoch, t 

o 

initial covariances, P. and P_ 

A B 

thrust/coast switching times, T , i=l,...,N 
final epoch, t 

r 

flag for covariance propagation method (state transition 
matrix or covariance integration) 
final state, X 

r 

final covariance, P A and P._ 

A B 

transition matrix, <f>(t_, T ) 

1 J 

dynamic noise, Q(t ) 

r 


Remarks : 

C0VP will propagate two sets of covariances which are usually 
the assumed (filter) estimation error covariance and the true 
(real-world) estimation error covariance. The latter covariance 
propagation is done only if generalized covariance analysis is 
desired. Two propagation methods can be chosen: state transition 

matrices (standard option) and integration of covariance matrix 
differential equations which are described in more detail in PROP 
and PATH, respectively. 
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C0VP--3 
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5.3.2 Program 


DATA 


Purpose: To read and edit user input data 

Input: (external) 

Output; listing and/or error report of user input 

Remarks ; 

DATA reads the first record of STM file to insure that the augmenta- 
tion parameters to Be used in error analysis are a subset of the parameters 
used at STM generation. 
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5.3.3 Subroutine DYNO 

Purpose: DYNO computes an effective process noise covariance due 

to time varying thrust errors 
Input : o thrust parameter uncertainties , q 

o correlation time, t 
o interval length, AT 

o basic 6X6 state transition matrix, $^+1 k 
o thrust transformation matrix to rotate thrust_parameters 
into cartesian coordinates , h 

o spacecraft state at start and end of interval, (X^, X ) 
Output: o effective process noise matrix, k 

Remarks : 


The process noise model assumed is a stationary Gauss-Markov 
process. Since the direct evaluation of a process noise covariance 
from this mode_l is time consuming, DYNO computes an analytic, approxima- 
tion to the actual process noise. The justification for this may be 
found in Appendix 9.3, The equations used are as follows 

« Vi,k ■ '<*'■'> [»\ + i + w Vk + i,iJ 

where 


(2) 

At t k+l “ \ 


(3) 

Q 

II 

^ ro l - 

> 

n 

V 

H 



j^O » At < t 


(4) 

R(At, t) = h. t At 


(5) 

\ = 3Y k /8h Ct k ), 

v^ = S/C velocity at t^ 

(6) 

p (t ) = cov [n 

q k *- 

M 
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DYNO - 2 


Since the effective process noise model is invalid over thrusting 

discontinuities, DYNO assumes that logic exterior to itself has 

adjusted propagation intervals to guarantee that a thrust on/off 

event does not occur in the calling interval, (t, , t ) 

k k+i 
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DYNO-4 



RETURN 






5.3.4 Subroutine 
Purpose : 

Input : 


Output ; 


FILTER 

To update augmented state knowledge covariances 
at a measurement event. 

o knowledge covariances before event, denoted 
by superscript ^ ^ 
o observation matrices 

o covariance of measurement white noise, R 
o logic control flags 

oo Kalman-Schmidt . weighted least squares, 
or user-supplied algorithm 
oo true or assumed covariance update 
o gain matrix if current update is for true 
covariances , K 

o updated knowledge covariances denoted by 
superscript ^ 
o gain matrix 


Remarks : 

As in subroutine, PROP, all equations below are written for a 
true covariance update. Wherever differences between true and assumed 
updates constitute more than simply dropping ignore parameter terms 
out of an equation, the difference is noted in the logic flow. Timing 
subscripts are not included here since the entire filtering operation is 
accomplished at a single time point. 

Using the linear measurement model described for TRAKM, Section 
5,3.22, results in the following equations for a covariance update. 
Defining first the measurement residual matrix, J 
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FILTERS 


1) J = H A + H B + HD + H E + H F + R 

X s u v w 


where 

2 ) 


rp m _ rp _ rp rp 

A = P H + C H+C H+C H + C H 

X XS S XU U XV V XW w 


rn fT rp m m 

3) B - P H + C H + C ~ H + C H +C^H 

S S XS X su u sv v sw w 

rp fp Hi rn rp rn rp 

4) D = C H +C H +UH +C H + C ~ H 

XU X SU S O U UV V uw w 

rp rn rT 1 rp rp rp rp rp 

5) E = C~H + C H + C H +VK +C~H 

XV X SV S UV U O V uw w 

rp rp rp rp rp rp rp rp rp 

6) F = C ~ H + C H +C H +C H +WH 

XW X SW S UW U VW V o w 

R is the measurement white noise covariance. If the update is to be 

for assumed covariances, one of the gain matrix subroutines, KSGAIN, 

WLSGAN, or USRGAN, is called to compute the state and solve-for gain 

matrices - K and K . True covariance updates use the gains previously 
X s 

computed by the FILTER pass which updated assumed covariances. 

For the Kalman-Schmidt filter, the updates proceed as follows, 
when K x and are the state and solve-for parameter gains, respectively. 

+ - T 

7) P = P - K A 

x 

8) C + =C“-KB T 

XS XS X 

9) c + =c~-kd t 

XU XU X 

10 ) c + = c “ - k e t 

xv xv x 


+ - T 

P = P -KB 
s s s 

+ - T 

C = C - K D 
su su s 

+ - T 

C = C - K E 
sv sv s 

c + = c - 

uv uv 
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FILTER- 3 


If, however, the update is for true covariances or assumed covariances 
for any aigurithm other than Kalman -Schmidt , several of the above 
equations change. While equations 9 ,10 ,,12 ,13 ,14 do not change , equations 
7,8, and 11 become, respectively 


15) 

4- 

r - ti 

T T 

P' = 

L p - K x A J 

- AK + K JK 

X XX 

16) 

+ 

r - t 

"j T 

c 

C -KB 

- AK + K J K 

xs 

L XS X 

J S X ! 

17) 

4* 

r - f] 

T 

P 

P -KB 

- BK + K J K 

s 

S X J 

s s s 


For true covariances the following equations are added 

,T 


18) 

19) 

20 ) 
21 ) 


C + = C - K F J 

XW XW X 


+ - T 

C = C - K F 
sw sw s 

C + = C 
uw uw 


C + = C 
vw vw 


Note that equations 15,16,17 are identical to 7,8,11 with additive terms. 
Therefore the standard procedure is to execute 7,8,11, and add the necessary 
terms from 15,16,17 if updating true covariances or using any algorithm 
other than Kalman-Schmidt , 
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Logic Flow 
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2 


Complete state and 
solve-for update 
[eqns. 15-17] 


FILTERr-4 


















5.3.5 Subroutine 


GPRINT 


Purpose; To print true covariances and their correlation coefficients, 

and dynamic noise covariance. 

Input; o true augmented state 

o true dynamic noise covariance 

Output ; (external) 

Remarks : 

GPRINT operates on true statistics in a manner analogous tp PRINT'S 
operation on assumed statistics. 
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5.3.6 Subroutine GUIDM 


Purpose : 


Input : 


Output : 


To compute guidance correction requirements and/or update 
the control error covariance 
o guidance initiation time, t^ 
o guidance cutoff time, t 

o guidance type: impulsive, low thrust, none (update 

control covariance only) 

o true estimation error covariance, P (t^) 

o variation matrix of targets WRT state at t , | 

c 

o sensitivity matrix of state at t WRT thrust controls, S 

o control covariance epoch, t 

o control covariance, P (t ) 

c o 

o transition matrix, ^(t^t^) 

o spacecraft acceleration (a T and a ) and mass (m T and m ) 

I c I c 

at t^ and t , respectively, and exhaust velocity, c 

o execution errors for impulsive guidance: proportionality 

(cr ) and two pointing angles (a and a.) 
r p p 

o control covariance epoch, t Q 

o control covariance, P (t ) 

c o 


Remarks : 

Five thrust controls are allowed: thrust proportionality, two 

pointing angles, guidance initiation time and guidance termination 

time. Selective weighting of the controls (W . . . ,W ) distributes 

1 n 

the control correction accordingly. Whenever the number of controls 
(either AV or low thrust) exceeds the number of targets, the guidance 
correction algorithm minimizes the weighted control correction. Ensemble 
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control corrections are pessimistically sized by manner the state 
control (actual-reference) covariance with the guidance matrix. 

A low thrust "vernier" guidance maneuver is performed between initia- 
tion and termination times of a primary guidance maneuver. The 
vernier removes state error accumulated since initiation of primary 
guidance or since the last vernier. Whereas the post-maneuver 
control covariance is normally set equal to the propagated knowledge 
at guidance termination , for primary guidance with subsequent 
vernier(s) it is important to set the post-maneuver control covariance 
equal to the knowledge at guidance initiation . 

Impulsive, or AV, -guidance computes an approximate mean AV by 
the Hof fman -Young formula using the AV covariance. 
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GUIDM-3 










GUIDM-4 



141 








GUIDM-5 








GUIDM-6 
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5.3.7 Subroutine KSGAIN 

Purpose: To compute gain matrix for Kalman-Schmidt filter. 

Input: o measurement residual matrix, J 

o cross .-covariance of state with measurement 
residual, A 

o cross-covariance of solve-for parameters with 
measurement residual, B 

Output; o gain matrix partitions for state, K , and 

solve-for parameters, K 

s 

Remarks : 

The equations coded are : 

K = A J" 1 
x 

K = B J _1 
s 
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5.3.8 

Program 

MEAS 


Purpose : 

To control measurement event processing 


Input : 

o current time 



o measurement type 


o generalized covariance flag 
Output: o updated augmented state covariance matrices 


Logic flow: 


ENTER 


*_ 

TRAXM 


Compute Observation 
matrices 


MENO 


Compute measurement noise 
covariance 
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MEAS-2 






5.3.9 Subroutine 
Purpose : 

Input : 

Output : 
Remarks : 
According to 
from the input array 


MENO 

To return the measurement white noise covariance 
corresponding to the current data type, 
o current measurement code 
o array of measurement variances 
o measurement white noise covariance, R 

measurement code, MENO loads the relevant variances 
into the current R matrix. 
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5.3.1Q Subroutine 


PATE 


Purpose : 


Input : 


Output ? 


To control estate propagation and computation of 
transition matrix 

• current time, t 

IV 

• integration end time, t^^ 

• covariance integration flag 

• transition matrix from t, to t, .. 

k k+1 

• state vector at t, , , 

k+1 

• augmented state covariance at t ^ 








5.3.11 Subroutine 


PRED 


Purpose : 
Input ; 


Output : 


To predict covariance values at spme future time. 

• time predicted to 

• current time 

• true and assumed knowledge covariances 

• true and assumed knowledge covariances at 
predicted time 


Logic flow: 
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5.3.12 Subroutine PRINT 


Purpose To output state vector . covariances and their correlation 

coefficients, and state transition matrices for assumed 
statistics 

Input: o current state 

o current time 

o augmented state covariance matrix 


o state transition matrix 
o assumed dynamic noise covariance 
Outputs (external) 

Remarks : 

PRINT transforms data into user-oriented output, computes correlation 
coefficients, and writes results on an output file. 
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5.3.13 


SUBROUTINE PROP 


Purpose : 


Input : 


To map covariance matrices at time t^ to time 

t^ + ^ using the state transition matrix method. 

o covariances at t, after all event calculations 
k 

at that time , denoted by subscript k and 
superscript (+) 

o state transition matrix partitions over current 
time interval, denoted by subscript (k+l,k) 
o thrust parameter uncertainty 


Output ; 


o flag indicating propagation of true or assumed 
covariances 

o Covariances at t k+ ^ before events, denoted by 
subscript k+1 and superscript (*-) 


Remarks : 

Propagation of the augmented state covariance proceeds as 


Cl) 

where 

( 2 ) 


^k+1 Vt-l.k ‘hT+l.k + Q k+1>k 


P = 


xs 

T 

xu 

T 

xv 

T 

xw 


s 

T 

su 

T 

sv 

T 


o 

T 

uv 

T 


(3) 


0 

0 

0 

0 


0 

I 

0 

0 


0 

0 

I 

0 


0 

0 

0 

I 
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(4) 


Q = 


0 0 0 0 


0 0 


0 


0 0 
0 


0 0 0 0 

0 0 0 0 0 

0 0 0 0 0 

Combining equations 1-4 yields the following equations, where 
transition matrix subscripts are ignored. 


(.5) P 


k+1 


T T 

+ + + 

$ P, + 0 C +0 C 
k xs xs, xu xu, 
k k 


+ 0 C + 
xw xw. 


T 

$ + C 


xs, , . xs 
k+1 


+ C 9+C Q T ~ 

xu xu xw 9 + Q 

k+1 k+1 xw k+l,k 

T T 

(6) C” = $ C + +0 P + + 0 C + ^ 0 C + 
XS. . .. xs. xs s, xu su. + xw sw, 

k+1 k k k k 


(7) C~ = $ C + +0 C + +0 U + 0 c + 

XS SU^ XU O XW UW, 


u k+i xu k 


(8) C" = $ C + +0 C + +0 C + +0 C + 

XV k+l Xv k XS sv k xu uv k xw w k 


(9) C = $ C + +0 C + +0 C + +0 W 

xw^ + ^ xw^ xs sw^ xu uw^ xw o 


(10) 

P 

S k+1 

■ P T 
s k 

(11) 

C~ 

S \+l 

= c + 

su k 

(12) 

c“ 

sv k+l 

= c + 
sv k 

(13) 

c~ 

sw k+l 

= c + 

sw k 
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(14) 

C 

c + 

UV k+l 

uv 

(15) 


- 4 “ 

C" 

C 

UW k+l 

uw 


(16) 


C w k+1 " C vw k 


k 

k 


Note that all of the above equations include ignore parameter 
information, which appears only in true covariance propagation of 
generalized covariance analysis. The calling sequence to PR0P indicates 
whether the current propagation is true or assumed covariances. For 
assumed covariances, all equations and parts of equations deriving 
from ignore (w) parameters are not processed. The following flow 
diagram does not show this in detail, so an additional diagram is 
shown as an example of this ignore parameter by-pass logic. 

All matrix multiplications, additions and subrtraction are 

performed by calls to matrix operations routines. In order to avoid 

programming complexity, calling sequences to these routines are 

always executed. The logic to prevent unnecessary operations - 

for example, attempting to compute C , when there are no solve-for 

XS k+I 

parameters ~ is included within the matrix routines themselves rather 
than in PR0P. 
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Logic flow; 



154 

















5.3.14 Subroutine 
Purpose: 

Input : 


Output : 


PTRAN 

To generate state transition matrix partitions for 
dynamic parameters by numerical differencing 
o spacecraft state at beginning of interval, x 
o spacecraft state at end of interval, x^^ 
o interval length, At 
o parameter list 

o perturbation magnitude for each parameter 
o parameter transition matrix 


Remarks : 

The augmented state transition matrix, $ , may be subdivided as 


$ 0 


xs 


XU 


0 0 


xw 


f k+l ,k 


I 

0 

0 

0 


0 

I 

0 

0 


0 

0 

I 

0 


k+l,k 


where the subscript (k+l,k) refers to the time interval (t^, and 


V*i . = 8 \+l /3S k> solve - for 

k+l,k 


0 = 3x , /3u dynamic consider 

XU« . « i tv * X 

k+l,k 


. = 3 W 3 v ignore - 

k+l ,k 

The zero entry in the top row corresponds to the state's independence 
of measurement consider parameters over time transitions. All sensitivities 
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computed in PTRAN are by numerical differencing. All parameter transition 
matrix elements are computed by PTRAN unless the variational method for 
thrust is selected by the user. 

Note that no mention within the PTRAN flow diagram is ever made of 
solve-for.. consider or ignore status for parameters. This is done for two 
reasons. First, when the state transition matrix (STM) file is created, 
no such reference is needed because all parameter sensitivities are 
generated at once. Parameter type specification is made at error analysis 
execution time and may change from run to run. Second, if PTRAN is ever 
used to generate transition matrices in-line with filtering operations, 
it may be exercised in either of two ways. The first would give PTRAN 
a parameter list including, in order, the solve-for, consider, and ignore 
parameters. The transition matrix would be returned and partitioned as 
necessary for the filtering operations. Or, separate calls to PTRAN 
could be made for each of the solve-f or-consider and ignore options with 
their distinct parameter list. The extra time necessary for multiple 
calls is more than saved by eliminating logic necessary to distinguish 
parameter type within PTRAN. 
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Logic flow; 


PTRAN- 3 


Zero out parameter 
transition matrix 

Initiate parameter 
list counter to zero 


( > 

J 


1 1 j 




Parameter list VYES_ 
Completed ? / 


RETURN 


Is current 
parameter dynamic 


Is current parameter 
thrust related ? 


Load relevant 
column from thrust 
■transition matrix 
into current para- 
neter transition 
natrix 


/ Is current parameter 
v planetary ephem. or grav. 


! Is S/C within 3 \ 

sphere of influence 

\ radii b$|yp leVant / 


Locate and purturb 
parameter nominal value ' 


Compute perturbed 
trajectory 


[ Compute state 
sensitivity vector 
by differencing 
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P IRAN -4 




7 

Restore parameter 
nominal value 

— 

1 


[Load sensitivity 
vector into current 
parameter transition 
matrix 




5 . 3.15 


Subroutine 


SCHED 


Purpose : 
Input : 


Output : 


To find the next scheduled event 
o Trajectory end time, T 
o current trajectory time, T^ 
o number of events , n 

o array of event start times , stop times , time 

between occurrences of this event , and event code 


(T . T , AT, C) 

start’ stop’ ’ 

o Time of next event , T 


o next event type 


Logic flow; 
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5.3.16 Subroutine 
Purpose : 
Input: 


Output : 
Logic flow: 


SETEVN 

To output trajectory- information 

• augmented state covariance 

• state and parameter transition matrices 

• spacecraft state 
(external) 
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5.3.17 


Subroutine 
Purpose : 

Input : 


Output : 
Remarks : 

Let G = 0 + 


STAPRL 

To compute the negatives of the partials of the 
spacecraft state WRT station locations 
o station locations, (R,0»0) 

0 Earth obliquity, e 

o Earth rotation rate, <d 
o universal time from epoch, T 
o current time , t 

o partials of state WRT station locations 


oj(t-T) 


3X 
s 

3R 

3X 

s 

39 

3X 
s 

30 

3Y 

£ 

3R 

3Y 
£ 

39 

3Y 

£ 

30 

3Z 

s 

3R 

3Z 

s 

39 

3Z 

s 

30 

3X 

s 

3R 


=-cos9 cosG 

= R sinO cosG 
= R cos9 sinG 

= C sine sin9 + cose cos9 sinG) 

•= R cose sin9 sinG - R sine cos9 
= -R cose cos9 cosG 

= sing cos9 sinG - cose sin9 
•= - (R sine sin9 sinG + R cose cos9) 
= R sine cos9 cosG 
= m cos9 sinG 
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3X 
s 

30 

3X 
s 

30 

3Y 

s 

3R 

3Y 

s 

30 

3Y 

g 

3 ^ 

3Z 

g 

3lT 

3Z 


— wR sin© cosG 


= wR cos0 cosG 


= - O3COS0 COS £ cosG 


= toR cose sin0 cosG 


= toR cose cos0 sinG 


= wsine cos0 cosG 


s 

30 

3 Z 
s 

30 '" 


= -o)R sine sin0 cosG 


= -ojR sine cos0 sinG 
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5.3.16 Subroutine; 


STMGEN 


Purpose! To create an STM file containing the integrated state 

and augmented state transition matrix for all 
events except prediction events. 

Input ; o event schedules 

o final time 

o start time 

o STM file generation data as described in Error 
Analysis Functional I/O Section 3.2 
Output: (External) 

Remarks : 

The layout of the STM file is: 

oo record 0: array of parameter numbers augmented 

to the state at STM file generation (used by 
DATA for error checking) 
oo records 1 - Nr 

word 1 - event time 

2 - event type 

3- 8 - spacecraft state vector 

9-44 - state transition matrix 

45-224 - parameter transition matrix 
At STM file generation, no distinction is made between solve-for 
and consider parameters. A single parameter transition matrix is computed. 
Parameters will be separated in subroutine STMKDR. 
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5.3.19 Subroutine 
Purpose : 

Input : 

Output 


Remarks : 


STMRDR 

To read the STM file and prepare information 
for the error analysis- module 
o event time 

o event type 

o integrated state vector 

o state transition matrix to event time 

o parameter transition matrix to event time 
o integrated spacecraft variables 
see subroutine STMGEN for STM file layout 
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Logic flow 


STMRDR -2 




Store parameter transition matrix columns into solve- 
for, dynamic consider, measurement consider or, ignore 
covariance matrices as directed by user input. 


RETURN 
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5.3.20 Subroutine 


TARPRL 


Purpose : 

Input : 
Output ; 
Remarks : 
since 


To compute the partials of the spacecraft state 
WRT orbital elements, 
o orbital elements, (a ,e ,i 
o partials of state WRT orbital elements 


tan — = 


1 + e 


1 - e 


tan 


E 

2 


M = E - e sinE 

true anomaly can be expressed as a function 
v = v (e ,M) 

The evaluation of the desired partials can now proceed. The results 
are summarized below, 

a. Partials with respect to a. 


3x _ x 
8a a 


lZ = Z 

3a a 

8z _ z 
3a a 


b. Partials with respect to e. 


3x _ xq 3v 

3e r r 3e 


3y _ yq . 8v 


3e 


r + r 3e 


cosf2sin(w+v) - sinficos (w+v) cosi 


sinftsin(w+v) + cosficos (co+v) cosi 


8z = ZSL + r 


3e 


3e 


cos(o>+v) sini 


where q = 


r 2^2 

72t r - a - ae Cl + sin v) 


aeCl - e 2 ) 
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where 


c. Partials with respect to i. 
= r siih2sin(oj-K>) sini 
- r cosfisin(oi+\;) sini 


3 x 
3 i 

- 

3i 


3z 


r = r sin(w-H>) cosi 


d. Partials with respect to/1. 

■ y 


3x 

3fi 


^ = x 

an 

3 z = o 

an 

e. Partials with respect tow. 

cosfisin(«+v) - sinncos (to+v) cosi 
innsin(w+v) + cosncos (w+v) cosi 


3x 
3 ai 


[- 

j^-- sii 


3oi 

3z , . . . 

-- = r cos (w+v) sun 
3w 

f. Partials with respect to M. 

cosnsin(w+v) - sinKcos(w+v) cosi 
sinHsin (oi+v) + cosncos (w +v) cosij 
cos(u+v) sini 


ae sm v 
s = ^ 


3x 

xs 

+ r 3V 

3M " 

r 

+ 3M 

3 y _ 

?s 

+ r — - 

3M 

r 

3M 

3 z 

zs 

3 v 

3M " 

r 

r 3M 


[' 


‘] 


[l - sf 
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5.3.21 Program: 


TEAM 


Purpose: To control the execution of the error analysis module. 

Input: see Error Analysis Functional I/O Section 3.2 

Output: see Error Analysis Functionla I/O Section 3.2 

Remarks : 

TEAM performs only control logic functions. All analytic functions 
of the error analysis module are performed by routines subordinate to program 
TEAM, 

Logic flow: see Macrologic, Error Analysis Section 4.1.2 
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5.3.22 Subroutine 


TRAKM 


Purpose : 


Input : 


Output : 
Remarks ; 
Data types 


To compute sensitivities of current measurement 
type to the state and all parameters, 
o measurement code 
o spacecraft state vector, x 
o parameter list 

o observation matrices, (H , H , H , H , H ) 

x s u v w 

available . 

o earth based tracking 


oo 

2-way range 


oo 

2-way doppler (range-rate) 

oo 

3 -way range 


oo 

3<-way doppler 


oo 

differenced 2-way and 

3-way range 

oo 

differenced 2r-way and 

3-way doppler 

oo 

azimuth and elevation 

angles 


o spacecraft based tracking 

oo star-planet/target body angles 

oo planet limb angles (apparent planet diameter) 
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All Earth-based data types are usable for near Earth missions, 
and can be taken from any tracking stations desired by the user. 
Interplanetary missions use all Earth-based data types except azimuth 
and elevation angles. These types must also be taken only from Deep 
Space Network (DSN) stations. Nominally stored in the program are the 
locations for DSN stations Goldstone, Madrid, and Canberra, but these 
locations may be changed or others added up to a maximum total of nine 
stations. Spacecraft-based tracking is restricted to interplanetary 
missions . 

Given the measurement model 

Z = h (x, s, u, v, w) 

assuming linearity for small deviations from nominal: 

6y = H 6x + H 6s + H 6u + H 6v + H 6w 

— V — Q — 11 — \T T.7 


8 ^ax 

where H = •- , H 

x s 


3h/ 


3s 

— , etc. 


The Earth-based data types are modeled using the following 
definitions (see Figure 1 for geometry) 

= S/C heliocentric position and velocity 
= Earth heliocentric position and velocity 
= Station 1 geocentric position and velocity 
Station 2 geocentric position and velocity 
= S/C position and velocity relative to stationl 

= S/C position and velocity relative to station 2 

= Unit vectors defining direction of S/C from 

stations 1 and 2 respectively 


-V 


--E’ -E 


- 1 ’ -1 
-2 ’ 4 


— l ’ £-1 


-2 ’ —2 


2 - 1 * —2 
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For two-way tracking the following model and sensitivities result 


P = \L\ = Ij* - Ig - rj 


P = £ • u 


1) 9p/ 


3x 


9p/ 3(.V r ft ) = (u^, z T ) 


9 P/ 


9p 


. 3 ^>4> 


3s 


Cr.^) 


2) a p / 


, - 9p. 

3s 3x 


3 <> x > ij) 


3s, 


3) 0 P/3x = (8|5 V * 3 P/ 8 r ) 

— — n — ti 


4) 9p/ aT . _ ;T 
3 ^h~ “ P 1 


3 %/3r, 


5) 3P ^ = % 

— n 


3p/ 


_ 3p 


3 s. 3 (r ^ , fj) 


3 *4’ ^ 

3s, 


6) 3p / 


3p . 3 4’4 } 


3s_ 3 jc 3s^ 


For use in (4) above 


7) 3u 


1 / 3 ^ [- 1 ^’ 


3x3 
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Both data types also have bias terms : 

P = |pJ + 

T 

P = £. • u + 

s > so/ 8 b t - 1 

Observation types including three-way data, whether as is, or 
in differencing, are also know as QVLBI Cquasi-very long baseline inter- 
ferometry) data types. Three way data types are modeled as the sum of 
the two way types plus a timing error term for tanging and a frequency 
bias term for range-rate. 

9) p 3 = + p 2 + cAt 

10 ) p 3 = p 1 + p 2 + c — 

where At is the timing error, c the speed of light, and A ^/f the frequency 
bias term which results from drift error between the frequency standards 
at the two separate tracking stations. The sensitivity partials for the 
three way data types are formed by adding the partials computed for 
each station individually. The c At and c f terms are treated either 
as biases or part of the white noise term. The differenced data types 
are modeled: 

11) Ap = p x - P 2 - cAt 

12) Ap = p x - p 2 - c Af /f 

The partials for the differenced data types are formed by differencing 
the individual partials, with the following exception. Since 
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13) 8A P/3r =3Ap /3r 
— ti — r 


[%-%] 


and u^and u ^are very nearly equal (as are and P 2 ) for interplanetary 
missions, we use the following substitutions: 

14) A£ = r_ 2 - £ 


15) Ap = [u + u 2 ] ^ Ar_ 


1 + Hi^ 


16) % " H.2 = [^~ ~ Ap — 2 ] ^ p ] 


Spacecraft elevation is computed from 


17) 



If elevation is negative, a note is made to that effect on the output 
file, but the error is not fatal. 

For azimuth and elevation angle partials, since no velocity dependence 
occur, we let 

21 = geocentric ecliptic S/C position 

x g = geocentric ecliptic station position 

u = unit vector in x direction 
-s — s 

w = unit vector orthogonal to x^ and geocentric ecliptic axis 

01 = S/C azimuth, measured positive from north toward east (see Fig. 2) 

0 = S/C el -ation 

p = S/C range vector from station 

u = unit vector in £_ direction 

x = projection of p into plane normal to x 
—a — s 

u = unit vector in direction of x 
—a —a 
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21 ) x = x + psinB u 
—a — s — s 


22) ^ = ?a/|xj , w = ~ — T - £ 


(x + y ) 
s 


[• V V °] 


23) sina = u w 


o/ \ 3a T 3u 

24) cosa /3x = w — a/3x 


. 3u 3x 

du . —a . —a 

25) -a/3x = ~ — 

—a — 


26 ) = sin3 u ~ +pcosB u 9 ^/3x 

3x -s 3x — s — 


o- 7 \ 3a,. T 3u , T 3w /r . 

27) cosa /3x = w — a/3x + u — /3x 

— s — — s —a — s 


~ d u 

oo\ 9u m —a 

28) — a/3x = ~ — 

— S OX 

—a 


3x 

—a 

3x 
— s 


29) a/3x = I + psinB 9 ~s/3x + sin$ u 9 p /3x + pcosB u ^/3x 

S jXj — s — s — s — s — s 


30) 9 --/3x = 

— s 


w i W 2 


1 - W, 


W 2 - 1 


- W;L w 2 


2 2 * 
(x s + y s ) 


0 0 c 

To complete the station location partials we have 


31) 3(o, P)/ - iiSbll 

3s x 

~ — s 


3x 
— s 

3s 


This is the same form used for the range and range-rate partials, where 
3 x 

— s/3s^ comes from subroutine STAPRL (5,3.16). 

For spacecraft based measurements, we use the following definitions: 
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x = spacecraft position vector 
Xp = planet position vector 
p = x - x = planet range vector 

- -p — 

d = vector of star direction cosines 
p^ = vector of target planet orbital elements 
Rp = target planet radius 
y = star -planet angle 

6 = apparent planet diameter measurement - angle subtended by 

planet disc at the spacecraft 

z^ = zero vector, 3x1 

Star-planet angle: 

T 

32) cosy = d • u 

~ -p 


3y/ 


3x 


33) 3y /3x = 


siny 

1 


d T /3x 

■ . T T 
E - ij cosy 


— psiny 


34) 3Y /3i = z T 


35) > Y /3 E =f J- 


Apparent planet diameter: 
36) sin^/2 = R p/ 


3x 

jn 

3E 


_ lx 

3x 


3x_ 

3 E 


37) 


36 

3x 


2 R u 


38) 35 /3x = 0 


* 3x p / 3 p generated by TARPRL (5.3.20) 
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39) 


36 /3p - | 5 - 

— dX 

~P 


3x 

E=R 

8 

P 


• 

3 

x 


3x 

_EE 

3 

E 


The sensitivities of all angle measurements to their respective 
biases are unity, as were the corresponding sensitivities for range and 
range^rate , 

Several pl& ces in the foregoing derivations , the partial derivative 

of a unit vector is needed with respect to its "parent" vector. Therefore 

the following notations are made. Define u as a unit vector in the 

direction of a 

u = a / | | 

— | a | 

then 



then 




3u/ 9c = - /3a_. 


Also, TRAKM is designed to receive a parameter list of all solve-for, 

dynamic and measurement consider, and ignore parameters. For deter mi ning 

observation sensitivities to parameters, one matrix, H is defined 

P 

initally to include all sensitivities which are then later rearragned 

column by column into the H , H , H , H matrices. This simplifies 

s u v w 

the logic necessary to compute the observation sensitivities initially. 
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Logic flow: 


ENTER 


TRAKM'H 


Earth or S/C based 
Tracking? 


Earth 



EPHEM 


Compute Earth 
ephemeris 
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r ^ 
6 


TRAKM-14 


Star-Planet 


Compute planet 
limb partials (eqn 37,38) 
Load into H v 


Compute star-planet 
partials (eqns 33,34) 
Load into H„ 



Target ephemeris 
in parameter list? 












TRAKM-15 
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5.3.23 Subroutine 


USRGAN 


Purpose: To compute filter gain by user supplied algorithm. 

Remarks ; 

User must supply his own FORTRAN subroutine and see that it is 
compatible with the calling sequence in subroutine FILTER (5.3.4) 
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5.3.24 Subroutine 
Purpose ; 


Input : 


WLSGAN 

To compute gain matrix partitions for sequential 
weighted least squares filter, and maintain 
the non-consider covariances necessary for that 
computation. 

o flag indicating gain matrix computation, 
or non-consider covariance propagation 


o for gain matrix computation; 

oo observation matrices, H , H 

x s 

oo measurement noise covariance, R 
o for non-consider covariance propagation; 

state and solve-for transition matrices 
Output; o gain matrix partitions for state (K^) and solve- 

for parameters OO 


Remarks : 

The sequential, or recursive weighted least squares (WLS) 
algorithm implemented here in equivalent to a batch WLS filter if 
there is no process noise. Since process noise is a significant part 
of low thrust analysis, the WLS filter must be used recursively, 
because it has no batch equivalent. The sequential WLS consider filter 
acknowledge consider parameters only for covariance analysis, and not 
for gain matrix caluclation. Therefore a set of "non-consider" covariances 
for the state and solve-for parameters must be maintained at all times. 

This set also represents the filter analysis as it would be in non- 
consider form. 

Each time the knowledge covariances are propagated - except for 
prediction - subroutine PROP (5.3.13) also calls WLSGAN to propagate the 
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WLSGAN-2 


WLS non-consider covariances. At a measurement event WLSGAN computes 

the gain matrix partitions K and K , and also updates the non-consider 

X s 

covariances. All covariances in the equations below are non-consider, 
not knowledge. 


Non-consider covariance propagation: 

+T 


1) 

P~ = 

$p + + 

0 C 
xs 

2) 

C ~ = 
xs 

$C + 
xs 

+ 0 

xs 

3) 

P~ = P 
s 

, + 
s 



T - T 
$ + C 0 

xs xs 


+ 


Gain matrix computations : 


4) 

A = 

- T T 

PH + C H 


X xs s 

5) 

B = 

T -T T 

PH + C H 


s s XS X 

6) 

J = 

H A + H B + R 


x s 

7) 


-1 

K 

= AJ 

X 

-1 

8) 

K 

= BJ 

s 


Non- 

-consider covariance updat 

9) 

P + 

- T 

= P - K A 1 


X 

10) 


C + =C'-KB 



xs xs X 


+ - T 

P = P -KB 
s s s 


11 ) 



Logic Flow 


WLSGAN-3 
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5.3.25 Subroutine 


XGUID 


Purpose : 
Input : 


Output : 


To control the execution of a guidance event. 

• current time (guidance epoch) 

• time of last control epoch 

• guidance cutoff time 

• target variation matrix flag 

• finite burn flag 

• updated control covariance 

• target variation matrix 












5.4 Simulation Mode 


5.4.1 Subroutine 
Purpose : 

Input ; 


Output : 
Remarks : 


CSAMP 

To sample eigenvalue of a covariance and rotate 

back into state space and form a sampled vector 

state covariance, P 

reference state, X 

flag for determining sample option 

dimension of state covariance, N 

array- of eigenvectors (R) and eigenvalues (V) 

sampled state vector, X 

s 

array of eigenvectors and eigenvalues, R and V 

Each eigenvalue is sampled assuming a normal distribution 

with zero mean using the function ENUM 
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5.4.2 Subroutine: 
Purpose : 
Input : 
Output : 


DATAS 

To read input and initialize trajectory simulation mode 
See TSIM Program Description (Section 2.3) 
o nominal spacecraft state at all maneuver times 
o guidance or variation matrices for all maneuvers 
o nominal target condi tons 

o error distributions of trajectory and mission parameters 
o all parameters in proper units and reference frames 
o random number initialization sequence 


Remarks : 

DATAS will prepare all data for subsequent Monte Carlo operation. User 
options will specify the degree of data preparation necessary, e.g. , whether 
target variation matrices are input or should be computed and whether priori 
error statistics are available from a previous run. Guidance and variation 
matrices are computed as in Sections 5.3.6 and 5.4.9 respectively. 
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LOGIC FLOW: 


DATAS-2 













5.4.3 Subroutine 
Purpose : 

Input : 


Output ; 


Remarks ; 


GUIDS 

To design the control correction necessary for a guidance 
event in the trajectory simulation mode, 
o maneuver epoch 

o estimated spacecraft state at epoch. X 

III 

o reference state at epoch, 
o target body and stopping condition 
o target values and tolerances 
o guidance law: linear or non-linear 
o guidance policy: impulsive or low thrust 
o allowable thrust controls (U) and weighting (W) 
o target variation matrix Gji) or sensitivity matrix (S) 
o maximum number of iterations 
0 guidance matrix (.for linear impulsive AV) 
o design control correction (Av or AU) 
o estimated target error before and after maneuver 
o target variation or sensitivity matrix 
o cycle termination flag 


Non-linear guidance applies a linear algorithm in iterative fashion. A 
successful maneuver design occurs when the corrected trajectory meets all 
target conditions within their tolerances. A near successful design occurs when 
the corrected trajectory comes "close u to meeting the target tolerances. "Close" 
may be defined as some scaler of the target tolerances. Should the non-linear 
design sequence exceed the maximum number of iterations and not come close 
to target tolerances, then the maneuver is deemed hopless to make, and the 
mission cycle is terminated. The target matrix G(j or S) is recomputed only if 
target error has increased since the last iteration. 
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Logic Flow 


GUIDS-2 


ENTER 



Linear 


Guidance law 





non-linear 


TRAJ 

Propagate X to target and 
compute target error AT 

’ 




T < tolerances 



YES 


< 


NO 


Max iterations 




NO 


YES/ AT less than 


Impulsive 


previous 

iterations 


Guidance policy 


low 

thrust 


'AT less than 
previous 
, iterations 


INO 


NO 


More controls 
than targets 


YES 


PINV 


T T -1 
n=ifi (# ) 


NO 


1 

MAT 

r 

IN 

n=f 1 

! 

1 


More controls 
than targets 


vYES 



HO 

MATIN 

n=s 1 



YES 


I 

TARMAT 


I 

SCOMP 

Compute target variation 


Compute weighted 

„ . , 3T 

matrix it= ~ 


sensitivity matrix 

3v 




S = (an) w 


PINV 


T T -1 

n=s (ss ) 


Compute, control 


T 

Compute control 

covariance 


correction 

Av = n. AT 


AU = n. AT 

X E = X E + AV 




U = U + AU 
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GUIDS-3 
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5.4,4 Subroutine 


NOISE 


Purpose: To compute thrust acceleration perturbations due to time- 

varying noise. 

Input: o standard deviation in thrust proportionality 

and two pointing angles , a_ 
o noise correlation times, x 
o noise flag (yes or no) 
o time interval, Ajt 
o present discrete thrust error, Aa 

Output: o new thrust error, Aa_ 

Remarks : 

The form of the thrust noise assumes acceleration error components to 
be independent of each other and to have Gauss-Markov properties. Use is made 
of the function RNUM to find Gaussian zero-mean, unit variance random numbers. 

Logic Flow: 
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5.4.5 Function 


RNUM 


Purpose: To generate a Gaussian zero mean random number. 

Input: standard deviation, a 

Output : random number , I 

Remarks : 

There exist many random number algorithms all of which perform equally 
well. The method used here requires a CDC system function (RANF) which generates 
a uniformly distributed random number between + 1, 



RAN] 

! 

F 

Compute two r 
r^ and 

andom numbers 
C 2 




x = O COS (277-^) \ 

- 2 LOG 10 < r 2> 
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5.4.6 Subroutine SC0MP 

Purpose: To compute the sensitivity matrix of target parameters 

WRT control parameters. 

Input: o epoch, t 

o nominal state, X 

o 

o nominal target values, T^ 
o target cutoff condition 
o nominal control parameters , U 
o control perturbations, AU 
Output: o sensitivity matrix, S 

Remarks : 

The sensitivity matrix is computed by numerical differencing techniques. 
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5,4.7 Subroutine 
Purpose : 

Input : 


Output : 


SETUP 

To transfer real-world or nominal or estimated values 
into working arrays. 

o ephemeris and gravitational constants 
o spacecraft constants 
o thrust control constants 
o ephemeris and gravitational constants 
o spacecraft constants 
o thrust control constants 


Remarks ? 

SETUP is used to store appropriate constants into arrays which are 
accessed by other routines. For example, prior to designing a maneuver in the 
simulation mode, SETUP is called to insert nominal mission values so that 
GUIDS will design the maneuver under the proper assumptions. 

Logic Flow: None 
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Subroutine 


STAT 


Purpose: To compute cumulative mean and covariance of error 

vector 

Input: o Vector of actual values, X 

’ a 

o vector of reference values x 

r 

o vector dimension, N 
o number of previous samples, m 

o mean (X^) and covariance (C ) of previous samples 
Output : o mean (X) and covariance (C) of total samples 








5.4.9 Subroutine 
Purpose ; 
Input ; 


Output ; 


TARMAT 

To compute target variation matrix and target parameters, 
o initial time 

o initial state and thrust controls 
o target cutoff condition 
o target parameters list 
o state perturbations, AX 
o target parameters , T 

£ 

o target variation matrix, \p 


Remarks ; 

The variation matrix is computed using the product of the state transition 
matrix (from initial to target time) and a target transformation matrix (at 
target time) . 
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Logic plow; 


ENTER 


TARMAT-2 














5.4.10 Subroutine 


TSIM 


Purpose : 
Input : 
Output : 
Remarks : 
See TSIM 


To control overall logic of the trajectory simulation mode, 
see TSIM Program Description (Section 3.3) 
o printout and punched cards 

Program Description (Section 3.3) and Macrologic (Section 4.1.3), 
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TSIM-3 


Compute actual control corrections 
by sampling error distributions 
and adding to design correction 


SETUP 

Reset trajectory constants to 
actual values. 


Any more maneuvers 


Is target before next 
maneuver 


Propagate actual state 
to target 



Propagate actual state to 
target but do not update 
state 


Compute cumulative statis- 
tics for this maneuver 


Any more maneuvers 
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6 . PROGRAMMING GUIDELINES 

To facilitate the conversion of FRACAS to computer systems other 
than Control Data and to meet core restrictions certain programming 
practices should be followed: 

o Data names ; labeled common names, and routine names will be 
restricted to six or fewer characters 
o F0RMAT statement literals will be defined as Hollerith fields 
o Alphanumeric constants will be six or fewer characters 
o Numeric literals will have eight or more significant digits 
to force double precision on IBM systems 
o Input and Output will be in READ and WRITE statements using 
logical files 5 and 6 (TAPE5, TAPE6) 
o Variably sized matrices will be treated as vectors to comply 
with matrix routine requirements (see Section 5.1.1) 
o Files will be defined with the minimum allowable buffer size 
o All large arrays will be located in blank common in order to 
maximize memory utilization since blank common overlays the 
area used for program loading information at load time (this is 
CDC peculiar) 

° A labled common, for example W0RK, should be available for local 
array usage to save memory 

° Variables defining maximum array size should be in data statements 

in each subroutine to facilitate any subsequent increase (or decrease) 
in large array requirements 
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Labeled commons 
commons could be ; 
o FLAGS 
o ENGINE 
o C0NTRL 
o EPHMRS 
o INTGRT 
o STATE 
o EVENT 
o AUGMNT 
o TIME 
o WORK 
o C0NST 


will contain related data items. Some typical 

- contains all logic control flags 

- contains all low thrust engine parameters 

- contains thrust controls 

- contains planetary ephemerides and other planetary data 

- contains all trajectory integration data 

- contains spacecraft and planetary state vectors 

contains all event information 

- contains all state parameter augmentation data 

- contains start time, final time, current time, etc. 

- contains working -storage memory for local usage 

- contains commonly used natural constants 
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7 . FUTURE OPTIONS 


The inclusion of future options in FRACAS is simplified as much as 
possible by two main features: modularity, which has been emphasized 
previously, and program coding standards (discussed in Section 6. Programming 
Guidlines) . Potential options or changes to the program fall into seven 
categories . 

o Dynamic model - both S/C and environmental characteristics 
c Integration algorithm 

o Transition matrix generation techniques 
o Additional state parameters 
o New data types 
o New filtering algorithms 
o Guidance algorithms 

Changes to the dynamic model are localized to subroutine TRAJ. Of 
course, a major rewrite of the dynamics - most likely for different thruster 
modeling - would mean extensive changes to TRAJ. However, the only other 
important, effect would possibly be to change the state parameter list, which 
is discussed below. 

A new integration algorithm would force replacement of INTEG and 
possibly a reformulation of the dynamics in TRAJ, but this is again a localized 
change . 

Transition matrix generation is currently dene in TRAJ and PTRAN , both 
of which are controlled by PATH. Since these already accommodate both 
integration of variational equations and numerical differencing, the only 
other possible technique is an approximate analytic technique. Its implement- 
ation would require a new subroutine and a modified calling sequence in PATH. 


211 



Additional state parameters impact transition matrix generation 
(TRAJ, PTRAN) and both propagation (TRAJ, PROP) and updating (TRAKM, FILTER) 
in filtering. If transition matrices for these parameters were to be done 
variationally, or if covariance integration were to be used, the relevant 
equations would be needed in TRAJ. However, the existing PTRAN could compute 
numerically differenced transition matrices. PROP also is unaffected. It 
needs only a set of covariances and a set of transition matrices - 
all of whose dimensions are already variable . Since TRAKM 
currently evaluates all observation sensitivities analytically, new 
equations would be needed. If sensitivities by numerical differencing were 
preferred, a differencing routine could also be added. FILTER, like PROP, 
would be unaffected. 

One mere important possible effect of additional state parameters, however, 
would be felt throughout the entire program. If additional parameters resulted 
in a required increase of the maximum dimensions of any array, such as a priori 
consider covariances, a program-wide dimension change would be required. This 
is the reason for the programming guidline which requires that any time the 
maximum dimensions, rather than current working dimension, of an array are 
used by the program, those dimensions must be defined by a variable passed 
through common, and not by a local constant. 

New data types, once modeled, require only additional coding in TRAKM. 

A new filtering algorithm, assuming it is linear, could use the existing 
FILTER subroutine. Only the gain matrix calculations would be affected. The 
only exception to this would be either a batch 0 r non-linear algorithm. The 
batch algorithm is of questionable validity for the low thrust problem 
because of process noise. A non-linear algorithm violates the basic linearity 

assumption of FRACAS and would, in all likelihood, require a completely new 

program. 
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Guidance algorithms are, again, local to the GUIDM secondary overlay 
and to GUIDS (in TSIM), and have no macro- impac t . 

A potential problem area common to all future options is the over- 
running of the 70,000 o word core restriction. If the overrun is minor, 

o 

localized maximum array dimension could be reduced to gain core. However, 
if the overrun is major, only two alternatives exist. Either the 70K 
requirement must be abandoned, or the TEAM primary overlay must be divided 
into two or more primary overlays, which is certain to increase execution 
time and could also reduce some program flexibility, particularly with 
respect to event types such as guidance and prediction. 
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9. APPENDICES 


The following three appendices contain technical analyses in support 
of FRACAS program design. Each appendix is self contained with its 
own references. The first Appendix (9.1) discusses Error Sources for 
near-Earth and interplanetary missions. The major error source is due 
to thruster performance Appendix 9,2 is a study of numerical accuracy 
of the covariance formulation. For a pre-flight error analysis program 
using CDC 6000 series computers, the covariance form is sufficiently 
accurate with no need for a square-root formulation. Appendix 9.3 
evaluates the advantages of two different covariance propagation methods: 
mapping with transition matrices and integrating covariance matrix 
differential equations. Transition matrix mapping is recommended for a 
pre-flight program because of its computational speed. 
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APPENDIX. 


9.1 Guidance and Navigation Error Sources 


INTRODUCTION 


A necessary part of any mission analysis, in particular 
guidance and navigation studies, is the identification of all pertinent 
error sources. The following survey seeks to summarize those error 
sources which apply to near Earth and interplanetary unmanned missions. 
The emphasis is on missions using continuous low thrust propulsion, but 
results can be used in ballistic missions since they are a subset of 
low thrust missions. 

NOMINAL ACCELERATIONS 


Quite often error sources are given as some percentage of a 
nominal value. It thus becomes necessary to understand the relative 
differences among the various forces acting on a spacecraft. Figures 
1 and 2 illustrate, for the interplanetary and near-Earth environment, 
the major accelerations affecting spacecraft motioir. The range of 
low thrust acceleration covers about .5 lb to .01 lb thrustors. The 
values for radiation pressure (and drag) assume large solar arrays 
(area/mass ~1 m^) 


It is observed that for low altitude Earth orbits, the low- 
thrust propulsion system does not overcome drag deaccelerations 
until about 400 Km altitude. Furthermore, the thrust levels for 
near-Earth missions are much lower relative to primary body accelera- 
tion than for interplanetary missions, which means that many revolu- 
tions about Earth would be required to raise an orbit from low Earth 
altitudes to geosynchronous. Since nuclear electric power decays 
exponentially over long times (years) and not as a function of helio- 
centric distance, it is quite possible to have thrust levels greater 
than solar gravity for outer planet missions. 

I. DYNAMIC (NON-THRUST RELATED) ERRORS 


Radiation pressure - Errors from 1 to 3% (la) of the nominal 
radiation pressure are due to (1) surface degradiation as the thermal 
environment changes, (2) inability to predict radiative/absorptive 
properties of all materials involved, and (3) changing effective area 
due to attitude motion with respect to the sun. 
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ACCELERATION (KM/SEC ) 
















Gravitational - Planetary ephemeris determination and prediction is a 
function of the quantity and quality of Earth based observability, 
using both optical and radar instruments. Accuracy is usually about 
2 to 10 p-rad (.4 to 2 arc sec) with radial (Earth line of sight) 
better than out of plane better than transverse (along velocity) 
determinations. For terrestial planets (Mercury to Mars), RSS 
0 ~100 Km, outer planets a ~ 1000 Km, comets and asteroids 10,000 
Km. Mjiss uncertainties vary widely from 2% of nominal u for Pluto 
to 10 % for the Sun. Earth mass uncertainties are 10 _it % while 

comet and asteroid uncertainties are generally large, 1 to 50% of 
their nominal p. 


Venting - Semi-random accelerations can be caused by outgassing from 
various scientific instruments, RTG's , propulsive valve leakage, or 
attitude control miscoupling. These accelerations generally vary 
from 10 to 10 Km/sec 2 . 

Earth atmosphere - The upper atmospheric density varies with the 
time of day, solar cycle, and a host of other phenomena. These 
variations along with changing effective spacecraft area result in 
drag uncertainties about 1 to 10% of the nominal drag force. 

Asphericity - Planetary figures are dominated by the second zonal 
harmonic Jy. For Earth J2=10~3 ) moon J2=2xl0~^, Jupiter J2=3xl0 -2 . 
Except for the Earth, whose J 2 accuracy is about .01%, oblateness 
uncertainties for the planets and moon are about 10%. 


Miscellaneous - Accelerations due to solar wind, micrometeorite 
impact, general relativity, etc., are usually ignored because their 
aggregate acceleration is less than 10~ 13 Km/sec 2 . 


Table 1 summarizes the nominal accelerations and typical uncertainties 
for near -Earth (h=20,000 Km or 3 R e ) and interplanetary (r=lA.U.) 
regions . The dominance of one error source over another is only 
weakly related to their respective nominal accelerations . 


SOURCE 

ACCELERATION (km/sec 2 ) 


Nominal 

lo Error 

Earth gravity 

lxlO -3 

lxlO -9 

Solar gravity 

6x10 8 

6xl0~ 13 


-6 

-8 

Thrust 

1x10 

3x10 

Earth J 2 

3xl0 _8 

3xl0 -12 

Lunar gravity (max) 

3xl0 -8 

lxlO -12 

Radiation pressure 

lxlO" 9 

2xl0 -11 

Venting 

10 -13 

10 -13 

Miscellaneous 

io -13 

10 -13 

Earth atmosphere 

M 

O 

1 

-P' 

O 

10 -41 


TABLE 1. SPACECRAFT ACCELERATION ERRORS 




II. DYNAMIC (HIGH-THRUST RELATED) ERRORS 


High thrust or impulsive or chemical propulsion is characterized 
by small Isp (^v/300 sec) and short burn times. These propulsion 
systems are used for midcourse and/or orbit insertion. Pointing 
errors are associated with establishing and maintaining an inertial 
orientation during the burn. For Mariner class midcourse engines, 
lo pointing is about 7 m-rad. Proportionality errors result 
from propulsion parameter uncertainties and variations during the burn 
with a ~1%. Resolution or quantization errors are associated with 
cut-off sensing using timers and/or integrating accelerometers with 
a~.01 m/sec. Generally, proportionality and pointing errors 
decrease as burn length increases. 

III. DYNAMIC (LOW-THRUST RELATED) ERRORS 


The most popular thrustor by far is the electrostatic Mercury 
ion bombardment engine. Discussion will be confined to this thrustor 
type although the general technique in obtaining effective thrust 
errors can be applied to any other thrustor type. Figure 3 illustrates 
the typical configuration. The power conditioner moderates any power 
fluctuations from the solar array (or any other power source). 

Error sources are broken down into (1) accelerating voltage errors 
caused by voltage regulation, neutralizer variations, and local 
potential changes, (2) beam current errors are caused by mixture 
uncertainties among singly and doubly ionized Mercury, main vaporizer 
controls, and beam signal control, (3) beam spreading (with net 
resultant loss of thrust) is caused by distortions in electric field 
shape and physical grid warpage due to initial placement or on-going 
thermal effects, (4) pointing errors are caused by control loops for 
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FIGURE 3. MERCURY THRUSTOR 





automatically gimbaling and translating the thrustor array; use of 
thrustors in attitude control mode reduces the thrust in the desired 
direction and introduces normal forces, and (5) failures can be 
continuous arc-outs (shorts between the grids due primarily to 
impurities) and outright thrustor or power conditioner failure 
which require corrective action, either ground or spacecraft initiated; 
the time lapse between detection and correction may be significant if 
the ground is in the control loop . 

Combining all of the engine errors into effective thrust 
errors permits a general input to guidance and navigation error 
analysis programs . Table 2 shows the contribution of each error 
source to the total effective error. If each error is assumed independent, 

T = n l + ^2 / ) 1^ (1+6) where M/e = mass to 

n l + 2n 2 ' cos© / >1 e b charge of Mercury 


PARAMETER 

CALIBRATION 
ACCURACY (%) 

STEADY STATE 
Iff (%) 

CORRELATION 

TIME 

AT/ 

BIAS 

T (%) 

TIME-VARYING 

I, , beam current 

.5 

1.5 

Weeks 

.5 

1.5 

V^, voltage 

.5 

1 . 

Weeks 

.25 

.5 

cos9 , divergence 

2. 

3. 

Weeks 

2. 

3. 

n^, single ion eff. 

1 . 

5. 

Days -Weeks 

.02-. 05 

.1-.2 

ri 2 > double ion eff. 

20. 

25. 

Days -Weeks 

.5-1.25 

.5-1. 

6 , fudge factor 

30. 

30. 

Days -Weeks 

.15 

.15 

Pointing 

— 

2 deg 

? 

— 

3.5 cross axis 


TABLE 2. MERCURY THRUSTOR ERRORS 


the net bias is about 2% (la) and the time-varying thrust error (process noise) 
is about 3% and 2 deg. with correlation time about a week. The principal 
engine errors are beam divergence and pointing. 

II . MEASUREMENT ERRORS 


The primary data types for near-Earth and interplanetary missions are 
Earth-based, in particular range and doppler. Besides errors associated 
directly with the measurement. Earth based measurements are affected by 
station location errors. These errors include not only physical location 
errors but many other processes whose effect is to perturb the spacecraft/ 
station signal such that the station location seems to be in error. These 
effective station location error processes include polar motion, Earth 
rotation rate (which affects timing by UT-ET conversion), charged particles 
in both space and ionosphere, tropospheric refraction, and instrument 
related errors: signal delay, oscillator instability and synchronization. 
Figure 4 illustrates the various improvements in calibrating out error 
processes associated with effective DSN station location longitude errors. 
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FIGURE 4. DSN STATION LOCATION ERROR IMPROVEMENT 

Much of the longitude error is associated with timing errors which are 
common to all stations, thus longitude error is often correlated between 
stations at about 0.9. Although DSN location errors are on the order of 
2 meters, SPADAT (Earth satellite tracking network) and MSFN (Manned Space 
Flight Network ) have location errors about 50 meters. The higher location 
error present for near Earth missions is tolerated because of tbe shorter 
spacecraft to Earth distances and stronger "observability" of the Earth- 
based data types: range, doppler, and angles. 


Current range and doppler (range-rate) uncertainties for the DSN are 
shown in Table 3. Typically, range measurements are wieghted at a fairly 
pessemistic level of about 50 meters. A summary of current error levels 
illustrated in Table 4 for various Earth— based measurement systems. 

VHF using ground transponders and lasers are used for near-Earth tracking, as 
well as landmark tracking and range/range-rate from navigation satellites. 
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OBSERVER RELATED 
ERRORS 

UNCERTAINTY (la) 

RANGE (m) 

RANGE-RATE (mm/ sec) 

DSN locations 

3 

.1 

Earth rotation 

2 

.02 

Pole motion 

.7 

.03 

Ionosphere 

— 

.4 

Space plasma 

10 

.8 

Troposphere 

5 

.3 

Station equipment 

— 

.4 

— 


TABLE 3. DSN RANGE, RANGE-RATE ERRORS 


TRACKING SYSTEM 
(one r measurement/min) 

NOISE (la) 

BIAS (la) 



Range 

(m) 

r 

(mm/sec) 

Angles 
(m rad) 

Range 

(m) 

r 

(mm/sec) 

Angles 
(m rad ) 

DSN 

50 

1 


0 

0 


msfn/spadat 

10 

.7 

.8 

20 

10 

1.6 

LASER 

1.2 


.5 

.15 



.5 

VHF 

15 

100 

— 

30-100 

50-200 

— 


TABLE 4. EARTH BASED TRACKING ERRORS 
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One additional Earth-based data type is the use of near-simultaneous 
differenced range and/or range-rate from two tracking stations. This data type, 
called Quasi-Very Long Baseline Interferometry (QVLBI) is proposed for future 
near-Earth and interplanetary missions. Expected DSN QVLBI range differencing 
noise is 1 to 10 meters and .1 to mm/sec for range-rate differencing. Table 
5 illustrates the expected improvement in DSN QVLBI range measurements. 

Effective QVLBI measurements require improvements in the current DSN system, 
primarily in clock synchronization for range and oscillator stability for range- 
rate . 



Present capability, m 


Projected capability, m 

Error 

Simul- 

taneous 


Present 

Upper value 

Lower 

value 

toured 

simul. 

tanoous 

configuration 

Simul- 

taneous 

Near 

simul- 

taneous 

Simul- 

taneous 

Near 

simul- 

taneous 

Charged 

particles 

1* b 

' 

Faraday rotation 

0.1 

0.5 



Troposphere 

1° * d 

1 

Constant model 

0.5 

0.5 



Signal arrival 
time/ground 
delay 

Clock tync 

10* 

1000* 

10* 

1 

Mariner Mars 
1971 plan- 
etary systems 
3*s 

10 

10 

1 

1 

1 

Clock rate 
ol 1 AU 

3' 

3 

*b standard 
~10’ u 

0.3 

0.3 



Transponder 
delay in- 
stability 

0.1 

1 

Mariner Mars 
1971 

0.1 

1 




Projected 

configuration 


S-X down link, 
1976 

Historical data 
improved mop- 
ping. 1973 


Star source VLBI, 
1976 

H standard, 

1973 


TABLE 5. DSN DIFFERENCED RANGE ERRORS 

A useful data type when near the target body is on-board optical 
data in the form of star/target body angles, target limb angles, 
and natural satellite (if any)/target body angles. The most efficient 
system makes use of the already present TV imaging instrument rather 
than a separate navigation device. Table 6 illustrates the optical 
accuracy for three systems ,a current system (Mariner 9 with a 500 mm 
full length and 1. 1x1.4 deg FOV) , a projected system for outer planet 
missions (TOPS), and the Apollo on-board sextant. 


SYSTEM 


la UNCERTAINTY (urad) 


Noise 

Bias 

Distortion (constant) 

Apollo 

50 

50 

0 

Mariner 9 

75 

25 

25 

TOPS 

25 

10 

10 


TABLE 6. OPTICAL NAVIGATION ACCURACY 
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One further instrument which is often mentioned in low thrust 
missions is the accelerometer. However, no accelerometers exist; _ 

7 rr — / 

which can accurately measure the low thrust acceleration CIO -’gto 10 g) 
over long time spans. Both bias and scale factors are temperature 
dependent, about 10 g per deg. Farenheit. The experimental spacecraft, 
SERT II, used a "minature electrostatically suspended accelerometer" 
with a purported error of about 1%. 

SUMMARY AND REFERENCES 


A typicai set of error sources for the early 1980 's is illustrated 
in Table 7 for a Mercury orbiter SEP mission. The levels assume 
significant improvements over current values in almost every area. 


ERROR SOURCE 1c VALUE 

Initial 
S/C State 

RSS Position 

25000 Km 

RSS Velocity 

25 m/s 

EP 

Thrust 

Bias 

Magnitude 

0.5% 

Direction 

0.5 deg 

Noise 

Magnitude 

2.0% 

Direction 

0.5 deg 

Correlation time 

1 day 

Radiation Pressure 

1.5% 

Mercury 

Ephemeris 

In-plane (ecliptic) position 

5 Km 

Out-of-plane position 

15 Km 

Gravitational constant 

0.4% 

DSN 

Station 

Location 

Radius 

1 in 

Longitude 

.5 m 

Doppler noise (per 1 min.) 

1 mm/s 

Range noise 

1 m 

QVLBI 

Doppler 

.1 mm/s 

Range 

— ... 

10 m 


TABLE 7. 1980 MERCURY ORBITER ERROR SOURCES 
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More detailed discussion of all the aformentigned error sources 
can be found in the references (Table 8) . 


ERROR SOURCE 

REFERENCES 


Radiation pressure 

1,2 


Gravitational 

1,2,3,5,17 


Venting 

1,2 

Dynamic 

Earth atmosphere 

6,16 j 


Asphericity 

’ I 


Miscellaneous dynamic 

1,2 j 


High thrust 

6,8,9,10 


Earth based (near-Earth) 

4,11,12,14 | 

Measurement 

Earth based (interplanetary) 

4,7,11 


Optical 

14,15 


Accelerometers 



10,13 

TABLE 8. ERROR SOURCE REFERENCES 




P . E . Hong A 

PEH/ac 

ATTACH. 
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APPENDIX 


9.2 Covariance Accuracy 

9.2.1 Preliminary Results for Task 3 of Low Thrust OP Contract 


The original intent of Task 3 was to arrive at an analytic 
approximation to covariance accuracy as a function of all matrix 
inputs to the filtering algorithm. Several obstacles arose in 
the analysis , so that a reasonable solution could not be obtained 
even for the simplest filter configuration. Consequently, a new 
approach is now being taken. The approach and its results will 
be described in a later memo. 

Computational error on a digital computer results from 
having to express each number in the finite word length of the 
machine. A machine which expresses each floating point number in 
t binary bits can store only a t-bit approximation to the numbers 
desired. We therefore have two sources of error — one from the 
initial t-bit approximation to the actual number, and one from the 
rounding or truncation when computational results must also be 
stored within t-bits. Since our purpose is to investigate 
covariance ill-conditioning, we will not dwell on such items as 
integrator accuracy, or the accuracy of transition or observation 
matrices. Since any covariance is theoretically positive semi- 
definite, and the mathematical operations of filtering retain 
this property, anytime a covariance becomes indefinite it must 
result from numerical problems. Since semi -definiteness is a 
theoretical necessity, physical subtleties such as the fact that 
transition and observation matrices are not exact are irrelevant 
to the conditioning problem. Thus, the only significant error 
source is the accumulated error resulting from the individual 
operation of addition, subtraction, multiplication and division. 

The error bound associated with any single arithmetic 
operation on two numbers is easily determined. If the operands 
are exact, the resultant is accurate to within ± 1 in its least 
significant bit. Thus, the absolute value of the resultant 
relative error is less than or equal to 2 -t for a t-bit machine. 
However, if the operands themselves are in error, this result is 
no longer true in general. In particular, the subtraction of 
nearly equal quantities (or, equivalently, the addition of 
oppositely signed nearly equal magnitude quantities) can produce 
unreasonably large relative errors (see Appendix A) . 
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Since each element in the product of an m x k matrix with a 
k x n matrix results from the sum of k distinct products, the 
potential for calculating terms with large relative error always 
exists for the aforementioned reasons. Thus, even for a single 
element in the product it is impossible to define an absolute 
upper bound on the relative error. However, even assuming reason- 
able error bounds, if the error in each matrix operation is assumed 
to attain its maximum value, the computed relative error bound 
grows to unrealistic levels (see Appendix B) . 

The reason this bound grows excessively is the underlying 
assumption that for the Euclidean norm, defined in Appendix B and 
denoted by |j li e , the only guaranteed bound for the norm of 
a product, |/AB// £ , is the product of the norms, )jA|l £ li Sil & 

The impact of this assumption is particularly evident in the multi- 
plication f PI 1 , where $ is the state transition matrix and 

P the state error covariance. For short time propagations, 

U § p ti T ll - is approximately equal to j) P He , 

yet ii <£ He for a particular sample Mercury orbiter mission is 
about 1.6 * 1CP. For the Euclidean norm, i! ft ll^ - j] fy T i)^ 
so we predict a bound on the norm of f p & T ^to be 2.6 * 10 8 
times the norm of P, which is ridiculously large. For a simple 
comparison, assume we would like to estimate the error bound 
only for the propagation J p i T , using the same technique 
as in Appendix B. The straightforward analysis gives a relative 
error bound of 6 * 10“^, meaning that a single propagation reduces 
accuracy from about fifteen digits to five or six (assuming CDC 
6000 series single precision) . If we make one breakdown of $ and 
P from their original dimensions of 6 x 6 into four sub-blocks 
each of dimension 3x3, the relative error bound drops to a 
more reasonable estimate of losing about one decimal digit. 

However, the problem now becomes intractable analytically, since 
this same type of breakdown would be required for all consider 
parameters. Their magnitudes vary widely, as do the sensitivities 
of the state and observations to them and a special case would 
have to be made for each to compute meaningful error bounds. 

One alternate approach would use a more optimistic error 
bound where we assume the error is small relative to the computed 
product, i.e. 

llP hii ll £ * 2.' r HflB He. 

This bound is suggested by Wilkinson (Ref. 2, p. 84) as acceptable 
in certain cases, but it can easily be exceeded and is too 
optimistic for the kind of assurances we need from this study. 

It would also be only semi-analytic in that it would require the 
norms of all intermediate matrices in the filtering operation, 
rather than depending on a few inputs as defined in Appendix B. 
These additional norms would have to be obtained explicitly since 
there are too many to evaluate analytically-. The final option 
would be to assess the problem using statistical error bounds 
which is again outside the scope of this task. 
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APPENDIX A 


Potential Rounding Errors in Addition of Inexact Numbers 

Two examples are given here of problems which can be encountered 
in a finite word length machine. They are designed to be illustrative. 
They may seem extreme - they certainly do not represent a normal 
situation — but they can occur. For these examples we will 
assume for simplicity that our machine carries three decimal 
bits and it truncates rather than rounds. The first example produces 
an answer in error by 100% when all the input quantities are exact. 

The second shows an infinite relative error when the input numbers 
are not exact, but have been truncated after previous operations. 

Example 1 

Compute 

(.601) x (.427) + (.348) x (-.731) 

The result of the first multiplication is .254627 and the 
second is - .254388. Thus, the correct answer is .239 x 10 - 3 , 

However, the machine will truncate each of those products to .254 
and -.254 respectively, and the computed answer will be zero, 
yielding a 100% relative error. 

Example 2 

Compute (.23125) x (-.32) + (.121875) x (1.92) 

The actual magnitude of each product is ,234, but each multiplier 
which is larger than three digits must be truncated in the machine 
before it can be operated on. The computed results will then be 

(.731) x (-.32) = - .23392 = r - .233 

(.121) x (1.92'' = .23232 = T .232 

_2 

With these two results the computed sum will be -.1 x 10 which 
has an infinite relative error compared to the zero expected from 
the calculation. Similar results can be demonstrated with 
rounding, though the problem is slightly less severe than with 
truncation. 

On a machine like the CDC 6000 series , where we have nearly 
15 decimal digits it may take a considerable number of operations 
before relative errors become as large as those of the examples 
but they can occur in some situations. Also, relative errors much 
smaller than 100% can completely destroy the validity of any 
results . 
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APPENDIX B 


Reproduced from 
best available copy. 

Assume we would like to predict the state covariance error 
bound for a single filtering step- which consists of propagation 
between measurement times and an update at the latter measurement 
time. The resultant error depends on the initial error in the 
covariance; on the norms of the initial covariance, transition 
and process. noise matrices, and the observation and measurement 
noise matrices; and on how each matrix enters the calculations. 

til til 

Propagation from the k to the (k + 1) time point proceeds 







> C V 


The measurement update consists of 


A- t- i 




Pi.y,jL ti ' [ H P H T i R 

/K Khi t * 1 *" Va. *+l j 


where 


p 


/Lf 




P 


V 


K -. 


r 


A.. H/ 
Sl- 


p = state covariance at time t. after processing 
l/ /. measurements up to and including time t / 


Q « process noise matrix accumulated over interval t, 
A+ 1 . — *. k 


^ = state transition matrix from time t^ to t^ + ^ 


to t 


k+1 


iJ = observation sensitivity matrix at time t, 

'“a -h k+1 


(in 


observation noise matrix 


t- ■= measurement gain matrix 
'Vt-i 

Rewriting the filtering equations for the error bound analysis 
so that we have a single operation per equation and dropping 
unnecessary subscripts we have: 


A i ■ * 

A 2 * A 1 * T 

A 3 = A 2 + Q = P k+l/k 

A 4 - H A 3 

T 

A c - A. H 1 
5 4 
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a 6 = A 5 + R 
A 7 ' V 1 

A 8-\ T4 7- K *-, 
A 9 = A 8 A 4 


10 


A - A = P 

3 9 k+l/k+1 


From Wilkinson (Ref. 1, p. 115), we have for the error in 
the matrix multiplication AB. 


/it//. - n 2.' r ' it H //. // 5/4 


where 


E s= error matrix; difference between computed product AB 
and the actual product of AB 

///?//„- matrix Euclidean norm defined by g = f 

t * I 

n = matrix dimension common to A and B 


4 

£1 


.'4, 


i /a ~ 

u 


”• t — t* 

t 1 is defined from 2 1 = 1.06 2 where t is the number 

of binary bits allotted by the computer for each floating 
point number's mantissa. 


Similarly, for addition we can derive from equation 6.16 
(Ref. 1, p. 115) that 

\\£ik * x'- c [ilniU - iisiul 

We now make the following definitions 

A ± = correct matrix at the i step 

A^ = computed matrix at i^ step 

y th 

E^ = - A^ = total cumulative error at i step 

*t c -» If it £ 

& - il He 





e 


and for Q, J , H, R, 
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// <?//* 



} etc. 


We assume that p k/k ls ln error 6y V 11,2 errM: ln P k+l/k+l 

is denoted E 1f ., and becomes the E for the next iteration. Inputs 
are 6’ f , , , /? c « , R ;°the P matrix dimensions are n x n; 

and the measurement is an m-vector. 


The a priori error upper bounds at each, step are as follows 
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where f (m) is a simple first or second order polynominal in m 
determined by the inversion method selected. 


The resultant upper bound on ^ requires eight pages of 
algebraic manipulation to derive and^one full page to write, but 
the key term is the multiplier of 2 1; 











fe 



The relative error at the end of the ten steps, is 


I i E 10 '/ 

-< 1 0 / • F° r one sample trajectory on the approach phase of a 

Mercury* orbiter we have the following numerical values for the 
relevant norms and dimensions: 


p t 

k+l/k+1 l 


E 


or 


n = 6 
m = 1 

(3 *= 3 =23* IQ 3 

■10 ' p J 1U 

P'f = 1.6 * 10 4 

^ 

I s 7 - 1« 5 

Since all entries in the computation are positive, the actual 

resultant upper bound can be no less than the bound computed by 
analyzing only the 2~t, terms . Knowing from above that /o 

gives ' 10 ' P 




tv 




[ C ■ (/. 6 * to' 1 )*' t- 

(/)* {L c )^'(fo s J- [(l.t’W’ 1 ) + 

£ (/. 6 A /0 i ) i ~ {2,3 * /o 3 ) ■+ .3]) 
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Ignoring smaller terms 

”f„ /f„ * ( i '*■) (/° £ )- ((■ i */C’‘’ J Y J. 

- (>-*' )' 6. r */»**) 

For the machines- of current interest, the CDC 6QQQ series, t = 48, so 

2 _t l - (1.Q61 C2" 48 ). 

- 3.76 * IQ -15 

which, gives a final relative error bound of 

<0 //V ^ S' * /o' ° 

Having computed an error bound - which is ten orders of magnitude 
greater than the resultant desired matrix, we know that this 
particular technique cannot provide meaningful answers. 


i*/*> 
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APPENDIX 


9.2.2 Conclusion of Covariance Versus Square Root Filter Formulaiton Study 


The conditioning effects of process noise on knowledge covariance. 

The presence of process noise in the orbit determination algorithm 
has an important effect on resultant knowledge covariance ill-conditioning. 
While that presence does prevent ill-conditioning in some cases, no ill- 
conditioning is observed for the CDC 6000 series computers even in the 
absence of process noise. Therefore, the full covariance form of the Kalman- 
Schmidt algorithm is recommended. The square root formulation is not worth 
its expense. 

The analytical evaluation filter accuracy is not a feasible approach, 
as has been discussed previously (Ref. 1). The current memo describes a 
new approach and its supporting computer program. Since the key to filter 
accuracy is computer word length, we shall investigate the filter sensitivity 
to word size by simulating machines of different word length. The standard 
internal format for floating point numbers on digital computers breaks the 
word into two parts, one for the most significant figures in the number's 
binary representations (mantissa) and the other for an exponent. This format 
is exactly analogous to standard decimal scientific notation. The CDC 6000 
computers have a 60 binary bit word length, with 48 bits used for the mantissa. 
This is the longest word of any current production scientific computer and 
was therefore selected as the reference length for the current study. As will 
be shown later, covariance ill-conditioning is not a problem with this word 
length, which is certainly a prerequisite for its use as a reference. 

The simulation program, called BANANA (Bit Allocations Necessary for 
Accurate Navigational Analysis), was constructed with few modifications 
to existing software. Word lengths shorter than the standard 48 bit mantissa 
are simulated by truncating bits as necessary at the least significant end 
of the work after each arithmetic operation. This is performed in F0RTRAN 
by masking expressions, which are available on the CDC. Masking expressions 
are expressions in which logical operations are executed on the operands bit 
by bit. For example, if we would like to evaluate the effect of a 24 bit 
mantissa - equivalent to IBM 360 single precision - we define a masking 
variable as 36 binary ones followed by 24 zeros. Then, after each arithmetic 
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operation, we perform the logical product of the resultant with the mask. 

The leading 36 ones in the mask preserve the 12 bits of the exponent, and 
the first 24 bits of the mantissa, while the trailing 24 zeros blank out the 
24 least significant bits of the mantissa. Since the masking variable is 
defined by input, any word length may be simulated by the change of a single 
card. 


This masking operation would be very difficult to insert in the program 
if all filtering and propagation operations were coded in line. However, 
these equations are coded in the program G0DSEP as calls to matrix operations 
routines. Therefore the masking expressions need only be added in these sub- 
routines to simulate the shorter word length operations. As was mentioned 
in the previous memo (Ref. 1), the purpose of this task is not to evaluate 
numerical errors in the integrator, or transition and observation matrix 
generators. Consequently, each BANANA run does not regenerate a trajectory 
with transition and observation matrices each time. For each trajectory 
study, a single G0DSEP run is made to generate and store on tape all transi- 
tion and observation matrices. The observation matrices are for all possible 
data types which could be exercised in a given study, so that BANANA has a 
flexible measurement schedule. When BANANA reads the transition and observa- 
tion matrices required for a given comparison run, it truncates all matrix 
elements, so initial accuracy of each matrix is consistent with the word 
length being simulated. 

Given an infinite word length with which to perform all calculations, 
the orbit determination (OD) results will be exact (again, assuming all 
input matrices to be exact). Even with the grossly pessimistic error bounds 
discussed previously (Ref. 1), it is possible to determine a finite word 
length which would allow the computed solution to approach the exact solution 
to any specified accuracy. However, this word length would be prohibitively 
large. In general, as word length decreases, we will see a gradual divergence 
of the computed from the actual solution. Since this divergence is only im- 
portant as it affects our physical interpretation of the OD results, the mea- 
sure of divergence must reflect our knowledge of the physical problem. The 
primary criteria selected, then, are the orientation and dimensions of the 
position and velocity uncertainty ellipsoids as represented by the state co- 
variance. In order to compare these quantities we compute the eigen values 
and eigen vectors of the position and velocity 3x3 sub-blocks of the state 
covariance. The relative orientations are determined by computing the angles 
between corresponding eigen vectors. The dimensions are compared as relative 
in standard deviations. Both comparisons are. made to the 48 bid word 
length OD results. A secondary comparison is made between gain matrices. 

For a scalar data type, the gain matrix is a vector. Considering the parti- 
tions corresponding to position and velocity, each represents a vector. 

The comparison made, then, is the angle between the reference and the current 
gain matrix for that run. 

Tests for numerical instability in a given computing problem are often 
made by comparing the results of double and single precision computations. 

A simpler but less meaningful comparison is available on the CDC 6000 series 
through the F0RTRAN compiler. The compiler offers an option for either 
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truncating or rounding the results of arithmetic operations to fit them into 
the 48 bit word. Rounding effectively provides an additional half bit com- 
pared to truncation, and differences can often be observed between computa- 
tional sets in which all inputs and operations are identical except for the 
differences of rounding and truncation. The masking operation previously 
described makes no attempt to evaluate rounding - all computational results 
are assumed to be truncated. This masking does, however, allow comparison 
of word lengths close to, but more than one half bit away from the nominal 
48. Since numerical problems result from the random, cumulative effect of 
losing information in the least significant bits of each number, these errors 
accumulate differently according to word length. But this difference is 
extremely small for any two word lengths, both of which are sufficiently 
large that neither suffers from significant accuracy loss. The comparison 
made here, was between 48 and 44 bit mantissae, and in all cases the 44 bit 
results were deemed sufficiently close to the 48 bit that both maintain 
acceptable accuracy levels. An example of sensitivity to numerical error 
was found in a comparison of two runs in a region of numerical ill-conditioning. 
Two runs with identical a priori and one bit difference in word length yield 
covariances after one day of tracking which differ by more than an order of 
magnitude on the diagonal. 

Results : 

The sample trajectory selected was the last 40 days prior to Mercury 
sphere of influence (SOI) encounter for a 1980 Mercury orbiter. Although 
different trajectories were not studied, data types were varied to evaluate 
the effects of observability, and the process noise level changed to evaluate 
its impact on OD conditioning. The primary indications resulting from this 
study are: 

(1) the greater the disparity in observability among state vector com- 
ponents, the worse the ill-conditioning problem, and 

(2) moderate low thrust process noise levels have significant stabilizing 
effect on both long and short term OD results. 

A summary of the runs made may be found in Table 1. For those runs with 
process noise, the nominal 1<T error levels assumed were 27. in thrust magnitude, 
0.5° in thrust pointing angles and a one day correlation time. 


Table 1: BANANA Results Summary 


— -Filter Inputs 
Word Length ‘ — - — . 

QVLBI With 
Process Noise 

R, R Only With 
Process Noise 

R, R Only 
No Process Noise 

44 bits 

SEE 

TEXT 


36 bits 

X 

X 

Wild fluctuations in 
eigenvalues and eigen- 
vectors 

30 bits 

X 

X 

Negative eigenvalue at 
1 day 

29 bits 

X 

X 


27 bits 

Negative eigen- 
value at 5 days 

Negative eigen- 
value at 5 days 

X 
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The first and most important comparison for each case is the 44 bit 
to 48 bit runs. The worst case for this comparison was with conventional 
tracking (no QVLBI) and with no process noise. With comparative information 
available every two days during the arc we find that all axes of both position 
and velocity uncertainty ellipsoid to remain well within .1°. Standard de- 
viations remained well within .01 % relative error with the exception of a 
short period of time immediately following the first ranging point. There 
an error in the smallest position eigenvalue did reach . 27o (.1 °L in standard 
deviation) for two days. This eigenvalue was also seven orders of magnitude 
smaller than the other two. 

The conditioning effect of process noise was best indicated by the 27 
bit run with conventional tracking. After the first ranging point at five 
days, the smallest position eigenvalue went negative and remained negative 
for nearly six days. However, at the end of the tracking arc - no automatic 
stopping procedure was built in for negative eigenvalues - the angular dif- 
ferences of the position and velocity uncertainty ellipsoids from reference, 
were all less than one degree. All standard deviation errors were 3.5 °L or 
less. Thus, even though propagation of physically meaningless covariances 
occurred for over five days, the process noise level was sufficiently high 
to wipe out all a priori information. In other words, the latter part of 
the tracking arc does nothing but maintain a balance between knowledge un- 
certainty increases from process noise and decreases from tracking, with 
little effect from earlier information. We also note that ill-conditioning 
came later with process noise than without. In both the 29 and 30 bit runs 
without process noise, negative eigenvalues were observed after one day of 
tracking, compared to the more than five days for a shorter (27 bit) word 
length with process noise. 
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APPENDIX 


9.3 PD0T vs. PHI 


INTRODUCTION AND SUMMARY 


In any linear error analysis program a major component is the 
propagation of state error covariances from one event to the next event. 
Two methods are generally used: integration of the covariance matrix 

differential equations (PDOT) and covariance mapping, with transition 
matrices (PHI). This study compares PDOT vs. PHI for low thrust trajec- 
tories from the viewpoints of both modeling accuracy and computational 
time. A key part of the evaluation is the process noise model which is 
especially critical for low thrust missions. Generally PDOT offers 
greater modeling flexibility and accuracy but at the cost of increased 
run time. It is recommended that for a pre-flight error analysis program, 
the PHI method and a semi-empirical noise model be used (along with certain 
operational guidelines) because it is 2 to 3 times faster than PDOT while 
retaining sufficient accuracy. 

NOISE PROCESS 


Given the nonlinear equations of motion 

x = x (x , u, ir) (1) 

where x is the spacecraft position and velocity, u_ are constant dynamic 
parameters, and n. are time-varying thrust parameters, these equations 
can be linearized about a reference trajectory such that 

fix = f <5x^ + g <5a + h 5q (2) 

where f,g, h represent . sensitivity partials (or transformation matrices) 
and 6x, Sx , 6u_ , 6n_ are errors in the respective dynamic parameters . Whereas 
eqn. 1 describes motion of the deterministic reference trajectory, eqn. 

2 describes the propagation of trajectory deviations resulting from 
dynamic and a priori uncertainties . A linear error analysis is concerned 
with the propagation of state errors through the uncertain dynamic environ- 
ment as affected by such events as measurement/state update and guidance 
(trajectory correction) . Of particular interest is the behavior of the 
ensemble trajectory error P, 

P (t) = E £ 6x(t)6x T (.t) ] 
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For low thrust missions the dominant error spurce hy far is thrust 
error (Reference 1 J , both bias and time-varying. Since a good OD filter 
can estimate biases fairly- accurately (Reference 2 ), the critical 
problem becomes the modeling of time-varying thrust error, 6q and associated 
h. Desirable features of the noise process are that 6;n (1) have a zero 
mean, (2) be stationary in a wide sense, that is, 6_n_ (tj) and 5q_ (t 2 ) are 
related only by the interval At = | t,-ti | , and ( 3 ) be time correlated 
such that the correlation between 6n(t^) and 5rj(t2) is inversely 
proportional to At. A convenient, yet simple, mathematical model which 
fulfills these characteristics is the Gauss/Markov process, which for 
simplicity is described in the one-dimensional case, 

6f|(t) = - — fin(t) + q (3) 

E [< 5 n(t)] = 0 


E 


[dn(t 1 )6n(t 2 )] = 



E [q(t)J = 0 

E [q( tl ) q ( 1 2 ) ] = f 6( V t 2 ) 

PDOT 

Propagation of P by numerical integration of the matrix differential 
equations is a straightforward application of eqn. 2, 

• T 

P = FP + PF* + Q 

where P is the augmented error covariance containing the normal spacecraft 
state, dynamic biases and time-varying thrust errors. 




0 0 0 


0 0 0 


0 0 o 

0 0 2pE [ 6n6n T ] 

This formulation has been incorporated into the low thrust error 
analysis program, G0DSEP, by Wayne Ratliff and exercised on a SEP 
Mercury approach trajectory. Figure 1 illustrates the growth of 
spacecraft position error for various values of noise. The nominal la 
thrust error (N) is 2% in proportionality and .5 degrees in pointing and 
a spherical a_ priori state uncertainty of 10 Km in each position component 
and .1 m/sec in velocity. The limiting value of the noise process as 
correlation time (t) approaches infinity is, of course, a bias. An 
important characteristic of both bias and noise is that a reduction, 
for example by an order of magnitude, of thrust error results in almost 
an exact corresponding reduction in state error growth. This is 
reassuring for the analyst who can then scale linearly the effects of 
error propagation corresponding to any given thrust error. It is also 
interesting to note that for the nominal error level, which corresponds 
to projected levels in the 1980's, the effect of noise resembles a scaled 
bias. This is illustrated graphically in Figure 2 and numerically in 
Table 1 (A through D) . Indeed, an empirical formula can be derived to 
estimate an "effective" bias for corresponding correlated noise, 

0 BIAS ” C * 30 + -° 5 T) °N0ISE (T in d3ys) (5) 

One further point observed in Table 1 is the correlation of the "considered" 
thrust error with the state. For correlation times about one day, these 
correlations are quite small which indicate the relative independence 
of process noise with respect to state. Of course, as correlation time 
increases the process noise looks more and more like a bias and the correla- 
tions approach significant values. 

An estimate of computer run time shows that each eigenvector event 
takes approximately .32 to .38 sec and each day of integration requires 
.15 sec when thrust noise is augmented to the basic state and .08 sec 
when bias is augmented (integration step size ** .1 day). These values 
are somewhat pessimistic because the PDOT formulation is not fully 
optimized which particularly affects eigenvector time. A typical 42 
day propagation with 3 eigenvector events takes 6.2 sec with biases and 
7.6 sec with noise. 
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FIGURE 1. 


THRUST ERROR PROPAGATION 
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FIGURE 2 „ THRUST NOISE vs „ BIAS 
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TABLE 1 . A . STATE ERRORS AND CORRELATIONS DUE TO THRUST BIASES 


srn dev 

X 

Y 

z 

VX 

VY 


9. 8 4 86765 2E* 8 3 

1 . 00 00 00 00 






l.C C 461282E+04 

-.58714297 

l.OCCOOGOO 





5 .7534 1UC3E+03 

-. 2995 1991 

.28553389 

l.OOOCCOOD 




1 . 0 9 2 922 4 7F- C 2 

. 9161 2376 

-.5 24 33054 

-.273751 24 

1. 0 090 0 CS 3 



: . i5er25i?t:-C2 

-. 5 2 71 2 6 65 

.97765 142 

.27782550 

- .50776702 

1 . 0000 00 03 


5.969J552 T E-G3 

-. 30583832 

.28371392 

.92778515 

- . 32576376 

. 2Ce4 3582 

1 . 005 O 3 DO 


-. 39C2 2912 

.25879867 

.15473349 

- . 63215493 

. 29946109 

.25425903 


-. 10941933 

-.200 15 219 

-.00263148 

- . 1034 312 5 

-, 33571167 

-.00 875973 


. 03244149 

-.01944589 

. 35196078 

.05320495 

-. 01948229 

.56619333 


TABLE l.B. 

STATE ERRORS AND 

CORRELATIONS DUE TO 

THRUST NOISE (t=10) 
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STD DC '! 

X 

Y 

Z 

VX 

VY 

V Z 

CO 

4. 2751444 9E *63 

1 . 0 0 00 0 0 oo 







4.67594<J49E*-C3 

= . 5 4 862167 

1.00000000 






2 8 5 2 4 7 8956E 4 fl 3 

“o 29164656 

. 28605824 

1. DOOOO0 00 





4.61184034E-G3 

o 7955 11 56 

- o 3 30 77 3 C 1 

-.21963966 

l.OOOCCQOO 




4.99703737E-03 

- 0 534Q 3867 

o960 20 777 

•285J276? 

- . 39075486 

1 . 0 0 00 coco 



2o45652446E°53 

“o 27867356 

o 24443 462 

.62268769 

- .31607338 

.26731667 

1. OOOQfiCOO 



01924655 

.00271218 

.00638302 

- . 2C824715 

.02078355 

.47626072 


IN- PMC 

0012 8416 

-.00773286 

-.00300732 

- . 8096 850 0 

— • 0 8L.. 6154 

-.00009406 


007-PLMJi 

. 00162190 

-.00822396 

.01441349 

. 01759503 

-. 00 17 0494 

. 17 1721 C9 


TABLE l.C. STATE ERRORS AND CORRELATION DUE TO THRUST NOISE (a=N, T=l) 


STO DtV 

X 

Y 

1 

V X 

VY 

vz 

4.5 7308 3 8 7E + J 2 

1 . 0 0 00 0 0 ;■)•') 






5 • 1 4 9 i 8 0 5 0E + 0 2 

-.4254 7987 

1 .0 0 0 0 00 0 0 





2. 92 76112 4E+02 

-.2328 94 50 

.23050252 

1 . 0 00 0 00 :s 0 




4 .6 7 0 05 633£ -0 4 

. 7 8 87 95 80 

-.26573676 

- .165 09474 

1.000 0 0000 



5.3 5 41 1358E-0 4 

-. 42397601 

.96380493 

.23919814 

-.33069923 

1 . 0 00 0 0 0 00 


2.525 45 32 4E -0 4 

-.24919426 

.22974747 

.80607646 

-.30052178 

. 2564 8664 

1.00000000 


-. 0 1 79933-. 

.0 024632 1 

.00550495 

-. 2056431 5 

. 01939383 

.07417920 


- • 0 G 12 DO 64 

-.0 0 7022.4 5 

- . 0 00 0 C631 

-. 00956242 

- . 07882556 

-.00039150 


. 00 15 1630 

-.00020341 

. 012430 70 

.0173750 1 

-.00153085 

.16703435 


-TABLE l.D. STATE ERRORS AND CONDITIONS DUE TO THRUST NOISE (a=.lN, T=l) 
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PHI 


Propagation of P by transition matrices is probably the most popular 
method. It relies upon transition or mapping matrices which can be generated 
by numerical differencing or By- a procedure similar to PDOT, that is, 
integrating variational equations . 


P(t) =<$> (t,t o ) P(t Q )$ (t,t Q ) + Q(t) 


( 6 ) 


where x contains only the basic spacecraft state and dynamic biases, 


HI) 


$ = 


fa e 


o i 


with $ = F<J> and $(t ,t ) = I 
o * o 


rJ 

Q = 


<Kt,s 1 )h(s 1 )E [«n(s 1 )6;n T (s 2 ^ T (s 2 )<f> T (.t ,s 2 )d S;L d 2 


It is apparent that this method should be much faster than PDOT if only 
because of the smaller dimensional state. However, an explicit assumption 
is that the thrust noise <5n_ and state error Sx are independent. As we 
have seen in the PDOT results this is not always true, particularly for 
long correlation times . A further drawback, is the need to evaluate the 
double integral for Q which would require substantial computer 
time unless reasonable approximations can be made . 

The program GODSEP propagates P by the PHI method. The state transi- 
tion matrix <|>, is obtained by integrating variational equations and the 
dynamic transition matrix 9, is obtained by numerical differencing. 

GODSEP also contains a Q approximation (see also Ref. 3), 


Q = R(At ,t) aH(t)+Kt,t o )HCt o H T Ct,t o ) 


where At - t-t 


H(t) = 


0 hE [ 6n(t)6n_(t) ] h^ 


( 7 ) 
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R(At ,x) = Hx At 


o<. 


= 2 for At>x 
= 0 for At<x 


It is important to note that Q is computed only from the last event to the 
current event and not at each integration step. This semi-empirical 
model for Q essentially translates thrust noise into an effective AV 
covariance and projects the covariance to the current time. For short 
intervals (At<x) Q resembles a bias. 


To determine the accuracy of the PHI method (and Q) a 20 day propagation 
is compared with the corresponding integrated values for the same Mercury 
approach trajectory used in PDOT (x=l day) . The overall propagation time 
is divided into smaller intervals (At) to examine their cumulative effect 
at the end of 20 days. Figure 3 shows the results for At=20, 2, and .1 
days. It is apparent that PHI propagation accuracy is good for large 
intervals, but breaks down for small intervals (At<x) . One characteristic 
which is difficult to observe in Figure 3, except for At=2, is the 
pessimistic estimate occurring early and an optimistic error later. 


The results of Figure 3 must be interpreted with respect to program 
usage. Guidance and prediction events generally require at least 10 day 
propagations which is comforting from an accuracy viewpoint. Measurement 
events fall in the other extreme of less than .1 day propagation. However, 
because each measurement alters the covariance, the cumulative effect of 
the combined measurement/propagation process must be considered. Figure 4 
examines the effect of taking 10 days of measurements starting at the end 
of a 20 day propagation (to build up the a priori covariance at the start 
of tracking) . The measurements represent a typical tracking schedule 
including range, range-rate, and differenced range and range-rate (QVLBI) . 

When estimation uncertainties are compared in Figure 4 it is seen that the 
behavior is similar between PHI and PDOT although the plateaus are different. 
Some of the differences may be attributed to the different a_ priori 
covariances at the start of tracking. The discontinuity in estimation error 
at 25 days is caused by a somewhat optimistic ranging point. Table 2 
summarizes the results after 10 days of tracking. The state vs. noise 
correlations remain small enough such that the knowledge error with 
PHI and Q is sufficiently close to that of PDOT. The approximate noise 
model of Eqn. 7 was arrived at semi-empirically and was found to be best 
overall. It obviously is far from perfect. Results of other models are 
displayed in Table 3 which compares the eigenvectors and eigenvalues at 
the end of tracking. Position error is more accurately predicted than velocity 
for all the tested models. 
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TABLE 2. A. PDOT STATE ERRORS AFTER 10 DAYS OF TRACKING 
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TABLE 2.B. PHI STATE ERRORS AFTER 10 DAYS OF TRACKING 





/KJ 

3 

(see Eqn 7) 


Velocity 

\Largest Eigenvector 
Eigenvalue (m/sec) \ Component 


R 

a 

X 

y 

z 

X 

y 

z 


(At) 2 

0 

.249 

.079 

,274 

.822 

.995 

.821 


(At) 2 

2 

.414 

.136 

.392 

.835 

.994 

.830 

GODSEP -> 

JgxAt 

0 

,471 

.174 

,451 

.842 

.990 

.837 


^xAt 

1 

.641 

.243 

.588 

.992 

.988 

.996 


PDOT 

.486 

.111 

.648 

.959 

,997 

.956 


TABLE 3. NOISE MODEL (Q) COMPARISON AFTER IQ DAYS OF TRACKING 


25 




An operational consideration in favor of the PHI method is the' 
availability of sensitivity matrices, <j> and 0, for output, which provide 
the analyst with a great deal of information on the trajectory and error 
processes. The cf> and 6 matrices between each event can also be stored on 
tape to facilitate later parametric error analyses (as opposed to PDOT 
which must store the F matrix at least once per integration step). 

As far as computational time is concerned, PHI requires about .27 
to .38 sec per eigenvector event and .04/. 12 sec per day of integration 
without/with thrust biases. A typical 42 day propagation with 3 eigenvector 
events takes 2.5 sec, and 8.0 with biases. 10 days of tracking (91 measure- 
ments) consumes about 30.9 sec printing every measurement and 13.5 sec 
printing every fifth measurement. This compares with 25.7 sec (printing 
every fifth measurement) for PDOT with noise. 

CONCLUSIONS AND RECOMMENDATIONS 


o For correlation times of the order of 1 day, the PHI method is 

about 300% faster than PDOT. Sufficient accuracy is retained because 
of the small state vs. thrust noise correlations; 

o for correlation times about 10 days, the PHI method is about 100% 
faster than PDOT if noise is simulated by effective bias (Eqn.5); 

o numerical differencing for the dynamic transition matrix requires 

about 50% more time than integrating the variational equations^ 

However, numerical differencing is straightforward to employ and 
does not require analytical partials (required in the F matrix) ; 

o a great deal of time is spent in print routines, particularly 
eigenvector and measurement , because the integration interval 
must be reduced to correspond to the current event and a consider- 
able amount of data manipulation is needed to display the information 
properly; 

o the need for high accuracy noise modeling in Q at small propagation 
intervals is diminished by measurement processing effects; 

o operational usage favors the PHI method. 

For the above reasons it is recommended that a pre-flight error analysis 
program propagate covariances by transition matrices, with Q and biases to 
simulate thrust noise. $ should be generated by integrating variational 
equations. As desirable options, PDOT and numerical differencing for 4> 
should be available. 
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